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The  effective  design  of  digital  control  systems  often  requires  the  availability  of 
highly  accurate  discrete  representations  of  continuous-time  systems.  The  discretization  is 
done  physically  via  a sample-and-hold  mechanism;  hence,  the  challenge  for  the  engineer 
is  to  produce  a good  discrete  model  of  the  sampled  system. 

Recent  successful  discretization  methods  are  based  on  high  order  Taylor  series 
expansions  or  lower  order  Taylor  series  expansions  combined  with  the  Carleman 
linearization  procedure.  In  this  work  we  introduce  a new  approach  based  on  a 
modification  of  the  Bellman-Richardson  extension  of  the  Carleman  linearization.  A key 
idea  is  to  approximate  truncated  terms  as  linear  combinations  of  kept  terms  with 
coefficients  that  are  polynomial  functions  of  the  manipulated  input.  These  coefficients 
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are  determined  by  solving  a numerical  optimization  problem.  The  resulting  linear  system 
is  then  integrated  over  the  sampling  period  to  yield  the  discrete-time  model  after  a 
straightforward  substitution  of  variables. 

Numerical  studies  are  conducted  to  evaluate  the  new  method  using  two  test  cases 
available  in  the  literature,  namely  a one-dimensional  model  of  an  isothermal  continuous 
stirred  tank  reactor  (CSTR)  and  a two-dimensional  model  of  a nonisothermal  CSTR.  The 
results  demonstrate  that  the  new  method  gives  excellent  discretizations  with  very  low 
order  expansions.  The  new  method  outperforms  the  Carleman  discretization  method;  in 
fact,  the  first  order  implementation  of  the  new  method  yields  discretizations  comparable 
to  those  from  a second  order  Carleman. 

In  addition,  a new  method  for  determining  suboptimal  feedback  control  laws  for 
nonlinear  dynamical  systems  is  developed.  The  new  approach  introduces  a structure  for 
the  control  law  similar  to  that  used  in  neural  networks.  The  nonlinear  dynamical  system 
and  the  performance  measure  with  constraints  on  the  inputs  are  then  linearized  using  the 
Caleman  method.  Thus,  the  dynamic  optimization  problem  is  transformed  to  an  algebraic 
expression  that  is  optimized  to  yield  the  parameters  of  the  control  law.  The  new  method 
is  applied  to  a one-dimensional  example  available  in  the  literature.  The  performance  of 
the  synthesized  feedback  controller  using  the  proposed  method  is  very  close  to  that  of  the 
optimal  control  law. 
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CHAPTER  1 
INTRODUCTION 


The  development  of  a discretization  method  for  nonlinear  continuous-time 
systems  is  motivated  by  the  computer-based  implementation  of  advanced  control 
systems.  Therefore,  discretization  schemes  that  lead  to  accurate  and  simple  discrete-time 
models  are  essential  tools  for  the  design  of  effective  digital  control  systems. 

The  lack  of  general  exact  discretization  methods  for  nonlinear  systems  forces  the 
control  designer  to  choose  one  of  two  alternative  approaches.  The  first  alternative  is  to 
design  a continuous-time  controller  based  on  the  continuous  model  and  implement  it 
digitally  through  a sample-and-hold  element.  This  approach  is  adequate  for  small 
sampling  periods.  The  second  alternative  is  to  design  a discrete  controller  based  on  an 
approximate  discrete  model.  Such  a model  can  be  obtained  using  classical  methods  such 
as  the  Euler,  Heun,  and  Runge-Kutta  techniques  (Chapra  and  Canale,  1988).  A 
disadvantage  of  these  methods  is  that  their  accuracy  rapidly  deteriorates  as  the  sampling 
period  increases.  This  disadvantage  is  overcome  in  a recently  developed  method  based 
on  the  Carleman  linearization  (Svoronos  et  al,  1994).  This  method  yields  a nonlinear 
discrete-time  model  by  integrating  a Carleman  approximation  of  the  original  continuous- 
time model.  Kazantzis  and  Kravaris  (1993)  showed  that  accuracies  comparable  to  those 
of  the  Carleman-linearization-based  method  can  be  achieved  by  using  Taylor  series 
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expansions  (with  respect  to  time)  of  high  order.  The  disadvantage  is  that  the  resulting 
discrete-time  models  are  rather  complicated. 

In  Chapter  2 of  this  dissertation,  a new  method  that  extends  and  improves  the 
Carleman-linearization-based  method  is  presented.  It  utilizes  a modification  of  the 
Bellman-Richardson  extension  of  the  Carleman  linearization  (Bellman  and  Richardson, 
1963).  This  linearization  approximates  a polynomic  system  with  a higher  order  linear 
system  by  introducing  an  additional  variable  for  each  higher  order  monomial.  The  basis 
of  the  extension  proposed  for  a Ath  order  approximation  is  to  express  the  neglected  terms 
of  order  k+ 1 as  linear  combinations  of  available  lower-order  terms.  The  coefficients  are 
determined  by  minimizing  a performance  measure  based  on  the  relative  difference 
between  two  estimates  of  neglected  high-order  monomials.  One  estimate  involves  linear 
combinations  of  available  lower-order  terms,  and  the  other  involves  products  of  first- 
order  terms  raised  to  the  appropriate  power.  Two  numerical  examples  are  used  to 
illustrate  the  performance  of  the  new  method. 

In  Chapter  3,  a new  approach  for  designing  a nearly  optimal  feedback  control  law 
is  developed.  The  three  well-known  approaches  to  optimal  control  are  the  Hamilton- 
Jacobi-Bellman  (HJB)  approach,  the  Calculus  of  Variations,  and  Potryagin’s  Maximum 
principle  (Kirk,  1970).  Each  of  these  three  methods,  under  appropriate  assumptions, 
gives  enough  information  to  determine  the  optimal  control  policy.  The  most  direct 
approach  to  the  determination  of  both  the  optimal  and  the  suboptimal  control  laws  in 
feedback  form  is  the  HJB  theory.  The  optimal  control  obtained  by  the  HJB  method  is 
usually  a time-varying  feedback  form.  In  addition,  it  is  impossible  in  most  cases  to 
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obtain  an  analytic  form  for  the  optimal  control  law.  Thus,  the  designer  is  faced  with  the 
determination  of  a suboptimal  control  using  approximate  methods. 

Several  approximate  methods  have  been  proposed.  One  of  the  first  approaches  to 
the  solution  of  nonlinear  problems  was  introduced  by  Al’brekht  (1961)  who  studied 
analytic  systems  and  discovered  the  optimal  control  as  a formal  power  series  by 
considering  Liapunov  and  Chetaev  methods  in  stability  theory.  An  analogous  problem 
had  been  considered  by  Lukes  (1969).  Pearson  (1962),  Burghurt  (1969),  and  Weber  and 
Lapidus  (1971)  introduced  a new  approach  to  the  solution  of  nonlinear  problems  by  using 
properties  of  the  well-known  linear  quadratic  problem  to  approximate  the  optimal  control 
of  the  nonlinear  system.  This  approach  involves  solving  the  linear  quadratic  problem  at 
different  values  of  the  state  variables.  Burghurt  was  able  to  precompute  the  optimal 
feedback  gain  by  expanding  it  in  Taylor  series  about  the  same  reference  values  of  the 
state  variables. 

A second  approach  was  suggested  by  Durbeck  (1965)  and  Garrard  et  al.  (1967). 
This  approach  is  called  parameter  optimization  and  is  based  on  generating  approximate 
solutions  to  the  HJB  equation  of  dynamic  programming  to  produce  suboptimal  feedback 
control  laws.  Generally  these  methods  require  a great  deal  of  computation  and  gives 
somewhat  worse  performance  than  the  approximation  to  the  linear  quadratic  problem 
(Garrard  et  al.,  1967;  Burghurt,  1969). 

Koepcke  and  Lapidus  (1961)  extended  the  idea  of  Liapunov  stability  to  generate 
near  optimal  control  policy.  This  technique  have  been  applied  to  the  feedback  control  of 
a tubular  reactor  (Parodis  and  Perlmutter,  1966). 
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In  recent  years  there  has  been  a growing  interest  in  implementing  neural  networks 
(NNs)  to  obtain  suboptimal  feedback  controllers.  Nguyen  and  Widrow  (1990)  and  Goh 
and  Edwards  (1994)  first  obtain  the  optimal  control  law  and  then  use  NNs  to  approximate 
it  by  a feedback  controller.  Recently  Goh  and  Edwards  (1997)  extended  their  work  by 
presenting  a design  method  that  uses  NNs  for  the  synthesis  of  dynamic  control  laws  for 
nonlinear  systems  where  the  state  of  the  system  is  not  measured  or  reconstructible.  NNs 
have  also  been  used  to  design  discrete-time  controllers  for  quadratic  performance 
measures  (Park  et  al.,  1996;  He  et  al.,  1996).  Following  a different  approach,  Yoshida 
and  Loparo  (1989)  used  a power  series  expansion  and  Carleman  linearization  to 
approximately  solve  the  optimal  quadratic  regulator  problem. 

The  method  presented  in  Chapter  3 combines  Carleman-based  linearization  and 
the  structure  of  nodes  with  activation  functions  used  in  artificial  neural  networks  to 
transformed  the  dynamic  problem  into  an  algebraic  form  that  can  be  optimized  to  yield 
the  feedback  control  law.  The  proposed  approach  is  based  on  the  proven  fact  that  the 
Neural  Networks  representation  can  accurately  predict  any  arbitrarily  nonlinear 
relationship  between  inputs  and  outputs.  This  will  lead  to  the  design  of  nearly  optimal 


feedback  control  law. 


CHAPTER  2 

NEW  DISCRETIZATION  METHOD 


Section  2.1  of  this  chapter  discusses  a proposed  modification  to  the  Carleman- 
based  linearization  method.  A method  for  designing  an  appropriate  performance  measure 
is  presented  in  Section  2.2.  The  structure  of  the  approximating  matrix  and  the  properties 
of  the  discrete-time  model  are  discussed  in  Sections  2.3  and  2.4,  respectively.  Finally,  two 
examples  to  evaluate  the  method  are  presented  in  Section  2.5. 


where  x(t)  is  an  ^-dimensional  state  vector,  ur  is  the  input  vector  that  remains  constant  in 
the  interval t e[tr,tr+I),  and  the  state  derivative  map  satisfies  f(0,0)=0.  Furthermore,  the 
input  ur  is  known  to  belong  to  an  allowable  range  [u™,  u„ux]  that  includes  the  zero-input 
vector.  Under  the  assumption  that  the  system  is  analytic  at  the  point  x = 0 , the  right-hand 
side  of  (2. 1)  has  a Taylor-series  expansion 


2.1  Modification  to  the  Carleman  Method 


Consider  the  nonlinear  continuous-time  system 


i(t)=f(x(t),ur)  t e[tr,tr+]) 


(2.1) 


oo 


[j] 


(2.2) 


whereA,'0(ur)-f(0,ur),  A'u(ut)  = ^/^  etc.,  and  where  x'ui  is  the yth  Kronecker 


product,  i.e.,  x'[2]  =x®x , and  if  x=[x,  x2]T  then  x't2)=[x2  x,x2  x2x,  x2]T. 
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Defining  the  vector  of  monomials  that  includes  the  first  Kronecker-product  vector 

xk'=[x'  x'[2]...x'[k]]T  (2.3) 

and  then  differentiating  xk  with  respect  to  time  yields  a higher-dimensional  equivalent  to 

(2.2): 

ik  = Sk(ur)xk  + Sk(ur)x'tk+1) +gk(ur)  + 0(x'tk+2))  (2.4) 

where 


a;, 
A 20 

a;2 

A21 

A'k, 

A2k.2 

A;k " 
A2,k-i 

Aj.k+i 

A 2k 

.A«  . 

, and  g’k  = 

1 

o ► 

o 

1 

s;  = 

0 

a;0 

•"  A 3>k-3 

A3,k-2 

. §'k  = 

0 

0 

0 

> 

rr 

o 

a;,  . 

withA'  = AJj  <S> In,_,  + In  ® A;_Itjfor  i = 2,  3,  k,  and  j = 0,  1, ...,  k-i+1,  and  where In is 


the  nxn  identity  matrix. 

The  elimination  of  redundant  states  in  equation  (2.4)  can  be  easily  accomplished  by 
replacing  vectors  xk  and  x'[k+11  by  lower  dimensional  vectors  xk  and  xfk+1\  respectively, 
that  lump  into  single  terms  all  components  of  x'k  which  are  repetitions  of  the  same 
monomial.  For  example,  in  the  two-dimensional  case,  x' = [x,  x2  xf  x,x2  x2x,  x2]T 
is  reduced  to  x2=  [x,  x2  x2  x,x2  x2]T  and  x'[2]=[xf  x,x2  x2x,  x2]T  is  reduced  to 
xI2]  = [x2  x,x2  x2]T . The  remaining  operators  in  (2.4)  must  then  be  modified  by  adding 
columns  that  correspond  to  permutations  of  the  same  monomial  to  one  column  and  by 
omitting  the  corresponding  rows.  This  does  not  modify  the  structure  of  the  matrices 
defined  in  (2.5)  but  reduces  the  size  of  the  submatrices  A'j  and  0 (except  forAJ,  and 
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A,'0).  Hence,  matrices  S',S'k,  and  vector  g^  are  respectively  replaced  by  Sk  , Sk, 
and  gk , of  reduced  dimensions,  and  (2.4)  becomes 

ik  =Sk(ur)xk  +Sk(ur)xfk+1]  +gk(ur) + 0(x[k+21)  (2.6) 

with  Sk  , Sk , and  gk  given  by  (2.5)  replacing  A'  by  the  reduced  A;j . 

The  Carleman  (1932)  linearization  method,  which  has  been  applied  to  nonlinear 
problems  (cf  Krener,  1974;  Svoronos  et  al.,  1980;  Tsiligiannis  and  Lyberatos,  1987; 
Kowalski  and  Steeb,  1991),  approximates  the  original  nonlinear  system  (2.1)  by  a linear 
system  that  is  obtained  by  truncating  all  terms  of  order  higher  than  k,  to  yield 

ik=Sk(ur)ik+gk(ur)  (2.7) 

where ik  approximates  the  state  vector  xk.  An  extension  to  this  approach  is  motivated 
from  a study  by  Bellman  and  Richardson  (1963),  who  use  a simple  one-dimensional 
example  to  propose  the  so-called  self-consistent  method.  Instead  of  neglecting  terms  of 
order  k+ 1,  the  idea  is  to  approximate  them  as  a linear  combination  of  terms  of  order  k or 
lower.  All  terms  of  order  higher  than  k+\  are  neglected.  The  coefficients  of  the  linear 
combinations  are  determined  by  minimizing  a mean-square  error.  The  Bellman- 
Richardson  example  concerns  only  a system  with  no  manipulated  input.  Furthermore,  the 
determination  of  the  unknown  coefficients  depends  on  the  initial  values. 

The  new  method  proposed  in  this  dissertation  follows  the  spirit  of  the  Bellman- 
Richardson  example,  but  develops  a framework  for  more  general  applicability.  The 
method  explicitly  includes  the  input  variable  and  attempts  to  capture  the  terms  of  order 
k+ 1 through  the  approximation 
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(2.8) 


where  the  matrix  H(y,ur)  has  elements  hfJ  which  in  turn  are  functions  of  a vector  of 

coefficients  y and  of  the  input  vector  ur.  Specific  recommendations  for  the  functional 
form  of  H(y,ur)  are  given  in  Section  2.3.  Using  the  expression  (2.8)  in  system  (2.6) 
yields  the  approximate  model 


where  Fk(y,ur)  = Sk(ur)  + Sk(ur)H(y,ur).  Note  that  for  the  special  case  where 
H(y,ur)  = 0 (or  equivalently,  for  x[k+1)  =0),  the  new  approximation  (2.9)  reduces  to  the 
standard  Carleman  approximation  (2.7). 

Since  the  control  variable  is  held  constant  during  the  sampling  period,  eq. 
(2.9)  can  be  integrated  analytically  to  yield  the  discrete  model 


A procedure  for  obtaining  matrix  <D  (y , u r ) and  vector  8 (y , u r ) can  be  found  in  Svoronos 
et  al.  (1994).  The  desired  estimate  of  the  state  x(tr+])  is  contained  in  the  top  n rows  of 
the  left-hand  side  of  (2.10).  The  next  rows  contain  a vector  of  estimates  of  the  second 
order  terms,  and  so  on.  Extracting  the  first  n rows  of  (2.10)  and  replacing  the  estimate  of 


each  product  (i.e.,  x,x2  ) with  the  product  of  the  estimates  (i.e.,  x,x2),  yields  the 
following  Ath-order  discretized  system: 


*k=Fk(Y>Ur)ik+gk(Ur) 


(2.9) 


ik(tr+l)  = ^(Y.Ur)ik(tr)+5(r,Ur) 


(2.10) 


A 


k 


i(tr+i)=£  ^ii(r,ur)i[lI*(tr)+5,(y,ur) 


i=l 


(2.11) 
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where  <E>u(y  ,ur)  is  the  submatrix  formed  from  the  top  n rows  of  <J>(y  ,ur)  and  from  the 
columns  that  correspond  to  xw,  8,(y,ur)  is  a vector  that  consists  of  the  top  n rows  of 
8(y,ur),  and  x[‘r  is  the  estimate  of  x[l!  obtained  by  employing  products  of  estimates  of 

first  order  terms.  For  example,  x(2)*  = [(x,)2  x,x2  (x2)2]  . 

2.2  Performance  Measure 

In  order  to  obtain  the  parameter  vector  y,  a specific  criterion  must  be  set  as  a 
measure  of  the  goodness  of  the  approximation.  Expression  (2.8)  implies  that 

x[k+,r  = H(y  ,ur)xk  (2.12) 

where  x*  is  obtained  by  replacing  in  xk  estimates  of  product  terms  by  products  of 
estimates  of  first  order  terms.  For  example,  for  the  two-dimensional  case  ( n = 2)  with 
truncated  order  k = 2 

A A A A A 

i2  = [x,  x2  (x2)  (x,x2)  (x2)]T 
and 

i;=[x,  x2  (x,)2  x,x2  (x2)2]T 

For  the  initial  value  ik(tr)  = 0,  (2.12)  is  trivially  satisfied.  If  the  discretization 
procedure  were  exact,  (2.12)  would  also  be  satisfied  at  the  end  of  the  sampling 
period,  tr+1 , with  xk  (tr+1 ) calculated  using  the  discrete-time  model  (2. 1 1).  We  judge  the 
goodness  of  the  approximation  by  the  difference  of  the  left  and  right  hand  sides  of  (2.12) 
at  time  tr+1 , with  zero  initial  conditions,  for  values  of  the  input  uniformly  spread  over  the 
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allowable  interval  [Umin,  iim**]  This  motivates  the  definition  of  a relative-error  vector 
e ( u ) with  components 


e,(u) = 


i = 1,2,  ...,N 


(2.13) 


where  the  subindex  / denotes  the  ith  component  of  the  vector,  and  N is  the  dimension  of 
vector  x[k+1)\ 

Consider  now  a sequence  of  inputs  {u,}"0  = {u0,u,,  where  u0  = umm  , 

um  = um«  > a°d  where  the  values  in  the  sequence  are  equally  spaced.  Then  the  proposed 
performance  measure  is 

J(Y)  = £!>?(”<>  C-14) 

t=0  1=1 


and  the  optimal  vector  of  parameters  is  found  by  solving  the  optimization  problem 

min  J(y)  (2. 15) 

T 

The  parameters  obtained  through  (2.15)  minimize  a relative-error  measure  as 
defined  in  (2.13).  Note  that  a key  idea  is  that  the  difference  between  the  two  estimates  is 
the  driving  force  provided  through  the  relative-error  e(  u ).  It  should  be  remarked  that  an 
absolute-error  measure  such  as  e(u)  = i[k+11*  - H(y,u)x*  cannot  be  used  as  a substitute 
for  (2.13).  The  reason  is  that  an  absolute  error  can  be  minimized  by  driving  ik(tr+1)to 
zero.  A trade  off  then  ensues  because  the  relative  error  (2. 13)  cannot  tolerate  a zero  input 
value  in  the  sequence  {u,)”0  since  this  would  make  the  denominator  identically  zero. 
This  difficulty  is  circumvented  by  replacing  the  zero  input  with  a small  nonzero  value. 
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2.3  Structure  of  the  Approximating  Matrix  Function 

The  elements  in  H(y,ur)  are  functions  of  the  coefficients  y and  the  input  ur . A 
polynomic  dependence  of  each  element  of  H(y,ur)  on  ur would  be  a simple  and 
attractive  choice.  However,  as  the  order  of  the  polynomial  in  ur  increases,  more 
parameters  must  be  identified,  which  leads  to  several  problems:  (i)  as  the  number  of 
parameters  increases,  local  minima  are  introduced  in  the  optimization  problem,  (ii)  too 
many  parameters  can  lead  to  an  optimum,  i.e.,  to  a y vector  value,  that  does  very  well  for 
the  inputs  included  in  J but  poorly  for  other  input  values  (overfitting),  and  (iii)  as  the  size 
of  the  optimization  problem  increases,  the  computational  cost  increases.  Therefore,  there 
is  a need  for  a structure  of  H(y  ,ur ) that  will  not  explode  the  number  of  parameters  when 
the  dimensionality  of  the  system  (n)  or  the  order  of  approximation  ( k ) increases. 

A first  step  is  to  set  the  trailing  columns  of  H(y,ur)  equal  to  zero,  starting  with 
column  rt+ 1 . This  is  equivalent  to  including  only  the  first-order  terms  in  xk  in  (2.8).  The 

following  reasoning  motivates  this  choice. 

An  examination  of  (2.4)  reveals  that,  after  the  truncation,  the  approximating 
dynamics  of  the  £th  order  terms  involve  only  three  nonzero  coefficient  matrices  (Ako,  Aki, 
and  Ato)  The  dynamics  of  the  k-\  order  terms  involve  four  nonzero  coefficient  matrices, 
and  the  number  of  nonzero  coefficient  matrices  increases  as  the  order  of  the  term 
decreases.  In  fact,  the  lower  the  order  of  the  term,  the  higher  is  the  order  of  the 
expansion.  This  suggests  that  the  first-order  terms  are  approximated  better  than  higher 


order  terms. 
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The  second  step  simplifies  H(y,ur)  further.  Consider  the yth  row  of  H(y,ur)  in 

(2.8) 

xf‘x2J...x^"  = hjjXj  + hj2x2+...+hjnxn  (2.16) 

where  we  have  used  the  fact  that  the  elements  hjn+1,hjn+2,...  have  been  set  to  zero  and 
where  the  integer  exponents  d,,i  = 1,2,...,  n,  are  given  by  the  yth  component  of  x[k+11 . If 
the  state  x, does  not  appear  in  xf‘x2J...x^  (i.e.,  d,  = 0),  it  does  not  directly  influence  the 
monomial.  We  therefore  set  the  corresponding  coefficient  in  the  right-hand  side  of  (2.16) 
equal  to  zero,  i.e., 

hjj  =0  if  dj  =0;  i = 1,2 #i . (2.17) 

Let  the  integer  t denote  the  number  of  coefficients  of  the  yth  row  that  are  set  to  zero  in 
this  step. 

The  last  step  provides  structure  to  the  nonzero  elements  in  H(y,ur).  Since  the 
performance  measure  involves  zero  initial  conditions,  x;(tr+i)  is  a function  of  ur  and  can 
therefore  be  approximated  by  an  coth  order  polynomic  expansion,  i.e., 

Xj=^(ur)  i = l,2,...,w  (2.18) 

where 


^i(“t)=rilu, + ri;u;j|+...+ri,u["1 


Then  an  approximation  of  the  monomial  that  is  first  order  in  xm  is 


Yd‘  Ydj  Yd”  ~J=d>Pd2  Pdm-1  e C dm+!  cd„ 

A]  A2  ...An  _(,]  (,2  •••Sm-1  Sm  Am  'sm+l  • • 


(2.19) 


(2.20) 


Such  an  approximation  can  be  written  for  each  of  the  (n  - £)  nonzero  terms  x,  that  appear 
in  the  right-hand  side  of  (2.16).  We  approximate  the  monomial  of  the  left-hand  side  of 
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(2.16)  by  the  average  of  all  such  approximations,  which  leads  to  the  following  final 
expression  for  the  nonzero  coefficients  of  H(y  ,ur ) : 


l __l_PdiPd2  tCi  fd.->|:d«i  Pd« 
njm  — , Si  S2  "'Wl  l»in  Sm+l",l»n  • 

1 n-i 


(2.21) 


For  example,  for  the  case  n = k = 2 and  where  ur  is  a scalar,  the  procedure  outlined  above 


yields  the  matrix 


H(y,u) 


0 0 0 0 
ooo 
0 o 0 

0 %\  0 0 0 


(2.22) 


where  t)l  = ynu  + y12u2+...  and  t,2  = y21u  + y22u2+.... 

This  structure  of  H(y,ur)  overcomes  the  problems  stated  earlier  and  yields  an 
excellent  approximation  as  shown  in  section  (2.4). 


2.4  The  Properties  of  the  Discrete-Time  Representation 
The  discrete-time  approximate  model  (2.11)  preserves  the  equilibrium  point  and 
the  local  asymptotic  stability  of  the  original  continuous-time  model  (2  .1). 

Theorem  1 . Preservation  of  equilibrium  point 

Let  (x°,u°)  = (0,0)  be  an  equilibrium  point  of  the  original  continuous-time  system 
(2.1).  Then  (x°,u°)is  also  an  equilibrium  point  for  the  discrete-time  approximate  system 
(2.11). 

Proof 
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Since  f(0,0)  = 0 then  A[0  = A10  = 0,  and  therefore  (2.5)  yields  gk  = 0 . From 
this  follows  that  8 = (y,0)and  8,  = (y,0)are  also  zero.  Hence,  it  directly  follows  from 
(2.1 1)  that  (0,0)  is  an  equilibrium  point  for  the  discrete-time  approximation. 

For  the  second  theorem  we  will  need  the  following  lemma: 

Lemma  1.  The  eigenvalues  of  A (reduced  system)  are  a subset  of  the  eigenvalues  of 

a;,. 

Proof 

For  the  £th  order  Carleman  linearization  the  reduced  matrix  Akl  can  be  related  to 
Akl  by  the  transformation 

A,,  = LA;,R  (2.23) 

where  L and  R are  elementary  matrices  with  the  property  that 

U„.R  = (2.24) 

and  q is  the  number  of  zero  columns  in  L.  For  example,  if  k = 2 and  n = 2 


0 

'l 

0 

O' 

'1 

0 

o' 

0 

1 

0 

0 

0 

1 

0 

0 

and  R = 

0 

1 

0 

0 

0 

1 

0 

0 

1 

Since  L and  R eliminate  redundancy,  the  column  of  L that  is  zero  corresponds  to  one  of 
the  two  rows  of  R that  are  identical.  For  k = 2,  R is  called  the  duplication  matrix  and  L 
the  elimination  matrix.  Explicit  expressions  for  constructing  them  are  given  in  Magnus 
(1988).  An  explicit  expression  for  R for  any  k and  n is 

R=SZ-S<eh  0e<2  ®-0e.k)eL.4 

>i=i  <i=i  *k=i 


(2.25) 
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where  e(| , e , ....  e,k  are  vectors  corresponding  to  the  /,,  /2,  ...,  ik  columns  of  the  In 
matrix  and  e , , is  a unit  vector  of  size  v(&,«)  - v(k  - 1, «) , where  (Svoronos  et  al.. 


1980) 


v(M)=  Z 


min  (t.n)/” 


i=l 


Furthermore,  e,^  Jk  is  the  same  unit  vector  for  all  permutations  of  i2,  /k , 


e.g., 


E112  — ei21  — E211 

For  the  form  of  eVj  Jk  that  satisfies  /, <i2 < ...</kthe  position  of  the  unity  element  is 
PV3..,k  i where 


ik-l  ik.j-l  n ifc.i— 1 n n 

^ X ^ -l^k  ^ ,1,2” Jk-3^k-2^k-l^k 

/k=l  /k.,=J  *k=l  ^k.2=l  *k-l=l  'k=l 

+•••+£ 

/,=!  iJ=l  lk=l 


(2.26) 


and 


6, 


if  /,  </2  </3  <-</k 

else 


(2.27) 


In  (2.26)  we  use  the  convention  that  ^*  = 0 if  3>a. 

P 

Once  R has  been  constructed,  L can  be  constructed  by  transposing  R and  changing 
any  repeated  column  to  zero  (i.e.,  if  / columns  of  RT  are  identical,  the  first  is  not  changed 
and  the  remaining  /-I  columns  are  replaced  by  the  zero  vector). 


From  (2.23)  and  (2.24)  it  follows  that 
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Akl-XInk_q=L(A',-XInk)R. 

Let 

Q = Ai,-Mnk 

then 

Aki  " = LQR  . 

If  A.  is  an  eigenvalue  of  Akl  and  w a corresponding  eigenvector 
(Ak,  -AInk_q)w  = 0 

hence 

LQv  = 0 (2.28) 

where  v = Rw. 

If  Qv  = 0,  then  A is  an  eigenvalue  of  Akl  and  v is  a corresponding  eigenvector. 
Given  the  structure  of  the  L matrix,  equation  (2.28)  implies  that  the  rows  of  Qv 
corresponding  to  nonzero  columns  of  L are  zero.  To  complete  the  proof  we  need  to  show 
that  the  remaining  rows  of  Qv  are  zero. 

The  definition  of  the  Kronecker  sum  of  a square  ^-dimensional  matrix  A and  a 
square  /w-dimensional  matrix  B is 
A®B  = A®Im  +I„  ®B 

and  since 

a;,  =a;,  ®i-i_  +i„®a;.u 

it  follows  that 


A2i  = An  © An 
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a;1=a;1©a;1©a;1 

A^,  = A[,  © Aj,  ©...©  Aj, 

>•  ' 

k terms 

also 

1.. 

1..  =ja,®i.®ij 

I k =-(In©I©...©In) 

n*  ^ \ p n “/ 

k terms 

therefore 

Q = A'kl -Ud,  = B©B©...©B 

where 

B = A'  --XI. 

11  k n 

Given  the  structure  of  R it  follows  by  direct  calculation  (see  appendix)  that  if  rows 
/ and  j of  R are  identical,  rows  / and  j of  [b©B©...©b]r  are  also  identical  for  any 
matrix  B.  For  example,  for  k = 2 and  n = 2 the  second  and  third  rows  of  R are  identical 
and  the  second  and  third  rows  of  [b©b]r  are  both  equal  to  [b21  (bu+ba)  b,J.  If 
rows  / and  j of  R identical,  so  will  be  the  / and  j rows  of  QRw  = Qv.  This  implies  that 
each  row  of  Qv  that  corresponds  to  a zero  column  of  L is  equal  to  a row  of  Qv  that 
corresponds  to  a nonzero  column  of  L (which  has  been  show  to  be  zero).  This  then 
implies  that  all  rows  of  Qv  are  zero,  and  therefore  \ is  an  eigenvalue  of  A', . 

Theorem  2.  Preservation  of  local  asymptotic  stability 
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If  the  original  continuous-time  model  is  locally  asymptotically  stable  at  the 
equilibrium  point  (x°, u°)  = (0,0) , so  is  the  discrete-time  approximation  (2. 11). 

Proof 

If  ur=  0 then  AJ0  = 0 and  A;o  = A;o  =... Ak0  = 0 from  the  formula  for  generating 
A'j . Therefore,  the  corresponding  reduced  matrices  A20,  A30,  ...Ak0  are  also  zeros.  It 
then  follows  that  Sk(0)  is  a block-triangular  matrix  with  submatrices  An,  A21,  ...Akl 
located  on  the  diagonal.  If  A.,  and  pjare  the  eigenvalues  of  AJj  and  A j , respectively, 
then  the  eigenvalues  of  their  Kronecker  sum  A'  are  given  by  A,  + Pj  (Graham,  1981). 
Similarly,  if  W;and  Vj  are  the  eigenvectors  of  A^and  A'_ljs  respectively,  then  the 
eigenvectors  of  A'  are  given  by  wj  ® v,  (cf.  Bellman,  1970  for  proof).  It  follows  from 
Tsiligiannis  and  Lyberatos  (1987)  that  the  eigenvalues  of  A',,  i = 2,  3,  ...  are  linear 
combinations  with  nonnegative  coefficients  of  the  eigenvalues  of  A(, . Since  the 
eigenvalues  of  Au  are  a subset  of  the  eigenvalues  of  A',  (Lemma  1)  it  follows  that  they 
also  are  linear  combinations  with  nonnegative  coefficients  of  the  eigenvalues  of  Aj,  which 
is  the  Jacobian  of  (2.1).  And  since  Sk(0)  is  block-triangular  with  An  on  the  diagonal,  it 
follows  that  all  eigenvalues  of  Sk(0)  have  negative  real  part  if  and  only  if  system  (2.1)  is 
locally  asymptotically  stable. 

Since  ur  = 0,  it  follows  from  (2.19)  and  (2.21)  that  H(y,0)  = 0,  which  further 
implies  that 

Fk(Y  ,0)  = Sk(0) 


and  from  (2.10) 
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<D(y  ,0)=  eFk(T,0)Tr . 

The  eigenvalues  of  Fk  (y,0)  determine  the  stability  of  the  discrete-time  model  since  the 
exponential  transformation  does  not  affect  the  stability  of  the  system.  Therefore,  if  the 
original  continuous-time  model  is  locally  asymptotically  stable  at  the  equilibrium  point 
(x°,u°)  = (0,0) , then  so  is  the  discrete-time  model. 

2.5  Examples 

The  performance  of  the  proposed  method  will  be  evaluated  by  employing  two 
previously  published  examples.  The  one-dimensional  example  is  first  presented  in 
Svoronos  et  al.  (1994)  and  then  in  Kazantzis  and  Kravaris  (1993).  The  two-dimensional 
example  is  presented  in  Kazantzis  and  Kravaris  (1993). 

(1)  One-dimensional  Example: 

Consider  an  isothermal  continuous  stirred  reactor  (CSTR)  with  the  prototype 
irreversible  reaction  2A— >B  with  rate  rA  = -2k,VcA,  so  that  the  system  is  described  by 
the  nonlinear  model 

V^=FcA_-FcA-2k,Vci  (2.29) 

where  V is  the  reactor  volume,  cA  is  the  concentration  of  species  A in  the  reactor,  cA  is 
the  feed  concentration  of  species  A,  F is  the  volumetric  flow  rate,  and  k,  is  the  reaction- 
rate  constant.  The  model  can  be  written  in  dimensionless  form  as  follows: 
dx 

— = -(l  + 2a)x  + au-ux-ax2  (2.30) 

dx 


where 


20 


X = 


u = 


F-F, 

F. 


2kVCA5 


a = 


= x. 


in 


and  where  the  subscript  ‘s’  denotes  the  reference  steady  state  variables.  Let  the  nominal 
steady  state  value  be  50%  conversion  which  corresponds  to  a=l. 

When  u = ur  (constant)  fort  e[xr,xr+1)  eq.  (2.30)  can  be  written  as 

— = -V,x-ax2  +V,  (2.31) 

dx 

where  V1  = l + 2a  + ur  and  V2  = aur.  The  exact  discrete-time  model  obtained  by 
integrating  (2.31)  analytically  is 


(2.32) 


where 


'2q  + V, +sV,T 


V2a  + Vj  - sj 


s = ylV ,2  +4aV2  . 


and 
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We  now  apply  the  proposed  method  by  developing  a first-  and  second-order 
realization  of  (2.6).  Using  a first-order  version  of  the  new  method  (k  = 1)  for  the 
linearization  of  (2.29)  leads  to  S,  = -V, , and  S,=ahn,  where  hn=£,  = 

ynur +Yi2ur +•••  • Hence  F,  = -(V, +a^,)  and  g,=V2  in  (2.9).  Therefore,  the 
corresponding  discrete-time  model  is 

x(T„, ) = <TL|Tx(t, ) + ^-(1  - e-L'T)  (2.33) 


where  Lj  = V,  +a^,. 

Using  a second-order  version  of  the  proposed  method  ( k = 2)  for  the  linearization 
of  (2.3 1)  leads  to  the  matrices 


S2  = 


-V,  -a 
2V2  -2V, 


and  S2  = 


0 0 
-2ahn  -2ah12 


where  hn  = £,2  and  h12  =0.  Hence 


f2 


-V,  -a 
_2(V2-a41)  -2Vj 


and 


g2  = 


resulting  in  the  discrete-time  model 
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x(xr+I)  = aV2 


(e  1 - 1)  (e  1 -1) 
a,  — : + a 


' 


+ a—. : t — -x2(xr) 


+ 


(Vi  +A.2)e  1 -(V,  +^,)e 
X2  -X, 


X,T 


x(Tr) 


(2.34) 


X2  — X, 


where 


a _ (V,+X2) 

a^-X,) 

t (V,+M 

2 a(X2  -X,) 


s2=V9V,2-4s, 


and 


Sj  = 2V,2  -2a2^,  + 2a V2 . 

In  what  follows  we  consider  two  cases  of  the  first-order  discretization: 
Case  1.1:  £,  =ynur 
Case  1.2:  =ynur  + yI2ur2  +y13ur3 


and  two  cases  of  the  second-order  version: 

Case  2.1:  =ynur 

Case 2.2:  =yuur  +y12ur2  + yI3ur3 
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To  determine  the  yij  values,  performance  measure  (2.14)  was  used  with  21  equally 
spaced  input  values:  {u/}™021  = {-0.9,- 0.8, 0.1, 0.000 1,0.1, 02,...,  1.2}  . This  results 
in  the  following  expression  for  £i  at  two  values  of  the  sampling  period  T : 

Case  1.1:  T=  0.05,  = 0.0462 ur 

T=  5.00,  =0.2901ur 

Case  1.2:  T=0.05,  £,=  0.0464 ur  -0.001 2 ur2  + 2.1583xl0-5  ur3 

T=  5.00,  =0.3405ur-0.1813ur2  +0.0760ur3 

Case  2.1:  T=  0.05,  ^ = 0.0462 ur 

T=  5.00,  =-0.2755ur 

Case  2.2:  T=0.05,  ^ = 0.0464 ur  - 0.0012 ur2  +2.0350xl0'5  ur3 

T=  5.00,  = 0.3397 ur  -0.1836ur2  +0.0792ur3 

To  test  the  discrete  model,  consider  the  input  sequence  ur  = 0.9(-l)r  for  xr  = rT , 
r = 0,  1, ....  9.  This  input  sequence  corresponds  to  alternatively  increasing  and  decreasing 
the  flow  rate  by  90%  of  the  reference  value,  which  represents  large  input  changes. 

Table  2-1  compares  the  output  of  the  exact  discretization  (2.32)  to  that  of  four 
other  methods,  namely  (i)  Euler,  (ii)  fourth-order  Runge-Kutta  (R-K),  (iii)  first-  and 
second-  order  Carleman  discretizations,  and  finally  (iv)  first-  and  second-order  versions  of 
the  proposed  method.  The  results  are  presented  for  sampling  periods  T = 0.05  and  T = 5 
(dimensionless  reference  residence  times).  The  table  shows  that  although  the  first-order 
version  of  the  new  method  (Case  1.1)  is  comparable  to  a first-order  Carleman  method  for 
the  lower  sampling  period  T = 0.05,  Case  1.1  is  more  accurate  for  the  larger  sampling 
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periods  T = 5 where  the  other  conventional  methods  (Euler  and  fourth-order  R-K)  fail 
completely.  Furthermore,  for  the  larger  sampling  period  Case  1 .2  (first-order  of  proposed 
method)  not  only  outperforms  the  first-order  Carleman  method  but  it  is  even  more 
accurate  than  the  second-order  Carleman  discretization.  Moreover,  both  cases  of  the 
second-order  of  the  proposed  method  prove  to  be  better  than  the  second-order  Carleman 
version. 

To  test  the  closed-loop  performance  of  controllers  designed  using  approximate 
discrete-time  models,  consider  the  deadbeat  response  to  the  following  triangularly  shaped 
set-point  sequence: 

f0.05r  for  r = l,2,...,5, 

M^)- jo  05(10  _r)  for  r = 6,7,...,10 

where  the  control  objective  is  to  reach  the  condition  x(xr+1)  = xjp(xr+1) . This  test  is 

particularly  relevant  because  an  exact  discretization  should  yield  a controller  that  meets 
the  control  objective  exactly.  Hence,  a good  discrete  approximation  should  produce  a 
closed-loop  response  that  follows  closely  the  triangular  set-point  signal.  Table  2-2  shows 
the  response  of  controllers  based  on  the  discrete-time  models  obtained  from  the  same  four 
methods  used  in  the  previous  open-loop  response  with  a sampling  period  T = 5.  The  R-K 
discretization  fails  to  yield  a solution  and  the  Euler-based  controller  performs  rather 
poorly.  The  controller  based  on  the  first -order  version  of  the  proposed  method  (Case  1.1) 
not  only  performs  better  than  that  of  the  first-order  Carleman  discretization,  but  it  is 
comparable  even  to  the  second-order  Carleman  version.  Note  also  that  Case  1.2 
considerably  improves  the  accuracy.  In  the  case  of  a controller  based  on  the  second-order 
of  the  new  method,  all  cases  yield  more  precise  results  than  those  of  the  second-order 
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Carleman  method,  especially  Case  2.2  where  the  closed-loop  system  follows  the  set -point 
extremely  well. 

(2)  Two-dimensional  Example: 

Consider  a nonisothermal  CSTR  with  the  reaction  A — > B,  with  rate 


rA  = -ZcAe  Re  where  Z is  the  frequency  factor,  E is  the  activation  energy,  R is  the 
universal  gas  constant,  0 is  the  reactor’s  temperature,  and  cA  is  the  concentration  of 


where  V is  the  reactor  volume,  F is  the  volumetric  flow  rate,  cA  is  the  concentration  of 
species  A in  the  feed,  (-AH)  is  the  heat  of  the  reaction,  p is  the  density  of  the  reacting 
mixture,  c is  the  heat  capacity  of  the  mixture,  0 _ is  the  temperature  of  the  feeding  stream 
and  Q is  the  rate  of  heat  input  to  the  reactor. 

The  dimensionless  form  of  (2.35)  can  be  written  as 


E 


species  A in  the  reactor.  The  system  is  described  by  the  nonlinear  continuous-time  model 


R0 


(2.35) 


— L = x„  -x,  -a(l  + Xi)e  ,+x’ 


b 


-x2 +ac(l  + x,)e  I+*2  +d(l  + u) 


(2.36) 


where 
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e-e. 


x2  = 


e. 


X 


li 


e.-e, 

e. 


u = 


Q-Q, 

Q. 


T 


ff 

V 


b 


c = 


d = 


(-AH) 

pces  A 

Q, 

FpcG, 


and  where  the  subscript  s denotes  the  steady  state  variables.  The  numerical  values  of  the 
parameters  used  in  this  study  are 

R = 8.314  kJ/kgmoleK 
E = 60,000  kJ / kgmole 
-AH=  10,000  kJ/kgmole, 

Z = 300,000  s'1 


p = 1 000  kg/  m 3 
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c = 4.0  kJ/kg  K 

and  the  operating  conditions  are  described  by 
cA^  = 10.0  kgmole/m3 

cao  = cas~ 4.286  kgmole/m3 
em=300  K 
0O=  0^=  400.08  K 
V = 0.01  m3 

and  F=  0.33x10^  m3/s. 

In  addition,  the  input  profile  is  the  same  as  that  of  the  one-dimensional  example 
considered  in  Section  3.1. 

The  Taylor  series  expansion  of  (2.36)  around  the  point  (xi,x2,u)  = (0,0,0)  is  as 
follows: 

x,  =-(l  + A)x,  -bAx2  - bAx,x2  +^-(2b  - b2)Ax2 

x2  = -cAx,  + (-1  + cbA)x2  + cbAx,x2  + ^c(b2  - 2b)Ax2  + du  (2.37) 

where  A = ae~b . 

The  first-order  approximation  (k  = 1)  using  the  proposed  method  yields 


F,  = 


-(l  + A)-bAh21 


cA  + bcAh21 


0 


and 
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where  h„  = h21  = Hz.  hzz  = Hi.  md  h33  = £2-  Moreover,  two  cases  of  the  first- 

order  version  are  considered: 

Case  1.1:^,=  y„ur  and  £2  = y21ur 


Case  1.2:  =yuur  +y12uj  and  £2  =y21ur  +y22u^ 

The  optimization  procedure  using  a sampling  period  T = 0.25  yields 
Case  1.1:  = -0.1235ur  and  £2  = 0.0489ur 

Case  1.2:  ^ = -0.1284ur -0.0408u?  and  £2  =0.0509ur  + 0.0012u* 


Accordingly,  the  discrete-time  model  of  (2.36)  is 
xi(Tr+i)  = T-^-r-[(an  -X,)ex’T  +(-a„  + X2)evr]x,(xr) 

A«2  A.  j 

+if^[eX,T-el,T]Xj(T') 


+ - 


‘12 


X2 


(e  a -1)  (e  1 -1) 
X2 


du 


*,(x„, ) = + eVr  h <*. > 

a,2(X2-Xi)  1 1 

+HrI<a"  'k,)el'T +('a" +x=>eX,,T]x.<T'> 

(a„ -X,Xel,T -0  (-all-X,Xew-l)' 


(2.38) 


X2-X^ 


du 


‘21 


where 


a„=-(l  + A)-ibA^;,  ajj^-bA-ibA^+itZb-b^A^, 

a21  = cA  + — bcA^2  , a22  = -l  + bcA  + — bcA^,  +— c(b2  -2b)A!;2, 
2 2 4 


and 
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_ (flu  + a22) V(^n  a22)  4(ana22  a21a12) 

1,2  - 2 

Using  a sampling  period  T = 0.25,  Table  2-3  compares  the  output  of  a fourth- 
order  R-K  simulation  of  (2.35)  to  that  of  the  Carleman  discretization  method  (first-  and 
second-order)  and  to  the  proposed  first-order  discretization  method.  The  results  show 
that  Case  1 . 1 of  the  proposed  technique  yields  considerably  better  approximation  than  the 
first-order  Carleman  method  and  is  comparable  (if  not  even  better)  to  the  second-order 
version  of  the  Carleman  technique.  Case  1.2  improves  the  accuracy,  though  only 
marginally. 
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Table  2-1  Open-loop  response  for  the  system  in  the  first  example 


Time 

Exact 

Euler 

R.K.4 

Carl.  1 

Carl.  2 

Case  1.1 

Case  1.2 

Case  2.1 

Case  2.2 

Sampling  period  = 

0.05  reference  residence  times 

1 

0.0409 

0.0450 

00409 

0.0409 

0.0409 

0.0408 

0.0408 

0.0409 

0.0409 

2 

-0.0060 

-0.0048 

-0.0053 

-0.0059 

-0.0060 

-0.0059 

-0.0059 

-0.0060 

-0.0060 

3 

0.0360 

0.0411 

0.0364 

0.0360 

0.0360 

0.0360 

0.0360 

0.0360 

0.0360 

4 

-0.0104 

-0.0083 

-0.0094 

-0.0103 

-0.0104 

-0.0103 

-0.0103 

-0.0104 

-0.0104 

5 

0.0324 

0.0383 

0.0328 

0.0324 

0.0324 

0.0324 

0.0324 

0.0324 

0.0324 

6 

-0.0136 

■0.0108 

-0.0126 

-0.0135 

-0.0136 

-0.0135 

-0.0135 

-0.0136 

-0.0136 

7 

0.0297 

0.0363 

0.0301 

0.0298 

0.0297 

0.0297 

0.0297 

0.0297 

0.0297 

8 

-0.0160 

-0.0126 

-0.0151 

-0.0159 

-0.0160 

-0.0159 

-0.0159 

-0.0160 

-0.0160 

9 

0.0277 

0.0349 

0.0279 

0.0278 

0.0277 

0.0278 

0.0278 

0.0277 

0.0277 

10 

-0.0178 

-0.0138 

-0.0171 

-0.0177 

-0.0178 

-0.0177 

-0.0177 

-0.0178 

-0.0178 

Sampling  period  = 

5.00  reference  residence  times 

1 

0.2185 

4.5000 

****** 

0.2308 

0.2179 

0.2163 

0.2187 

0.2187 

0.2185 

2 

-0.5952 

-148.5000 

****** 

-0.4286 

-0.5380 

-0.4894 

-0.5653 

-0.5475 

-0.5809 

3 

0.2185 

****** 

****** 

0.2308 

0.2179 

0.2163 

0.2187 

0.2187 

0.2185 

4 

-0.5952 

****** 

****** 

-0.4286 

-0.5380 

-0.4894 

-0.5653 

-0.5475 

-0.5809 

5 

0.2185 

****** 

****** 

0.2308 

0.2179 

0.2163 

0.2187 

0.2187 

0.2185 

6 

-0.5952 

****** 

****** 

-0.4286 

-0.5380 

-0.4894 

-0.5653 

-0.5475 

-0.5809 

7 

0.2185 

****** 

****** 

0.2308 

0.2179 

0.2163 

0.2187 

0.2187 

0.2185 

8 

-0.5952 

****** 

****** 

-0.4286 

-0.5380 

-0.4894 

-0.5653 

-0.5475 

-0.5809 

9 

0.2185 

****** 

****** 

0.2308 

0.2179 

0.2163 

0.2187 

0.2187 

0.2185 

10 

-0.5952 

****** 

****** 

-0.4286 

-0.5380 

-0.4894 

-0.5653 

-0.5475 

-0.5809 
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Table  2-2  Closed-loop  response  for  the  system  in  the  first  example 


Tim  Exact 


Sampling  period 

1 0.0500 

2 0.1000 

3 0.1500 

4 0.2000 

5 0.2500 

6 0.2000 

7 0.1500 

8 0.1000 

9 0.0500 

0.0000 


Euler 

R.K.4 

Carl.  1 

Carl.  2 

Case  1.1 

Case  1.2 

Case  2.1 

Case  2.2 

= 5.00  reference  residence  times 

0.0033 

***** 

0.0492 

0.0050 

0.0500 

0.0500 

0.0500 

0.0500 

0.0097 

***** 

0.0972 

0.1001 

0.1000 

0.1000 

0.1000 

0.1000 

0.0188 

***** 

0.1441 

0.1503 

0.1505 

0.1499 

0.1500 

0.1500 

0.0330 

***** 

0.1903 

0.2005 

0.2016 

0.1998 

0.1999 

0.2000 

0.0440 

***** 

0.2361 

0.2509 

0.2539 

0.2502 

0.24% 

0.2500 

0.0535 

***** 

0.1903 

0.2005 

0.2016 

0.1998 

0.1999 

0.2000 

0.0593 

***** 

0.1441 

0.1503 

0.1505 

0.1499 

0.1500 

0.1500 

0.0618 

***** 

0.0972 

0.1001 

0.1000 

0.1000 

0.1000 

0.1000 

0.0611 

***** 

0.0492 

0.0500 

0.0500 

0.0500 

0.0500 

0.0500 

0.0574 

***** 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 
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Table  2-3  Open-loop  response  for  the  system  in  the  second  example 
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CHAPTER  3 

NONLINEAR  OPTIMAL  CONTROL 


3.1  Statement  of  the  Problem 

Consider  the  nonlinear  optimal  control  problem  given  by  the  functional 
minimization 


lr 

min  J = M(x(tf ))+  J N(x(t),u(t))dt  (3.1) 

l0 

subject  to  the  dynamic  constraint 

x(t)  = f(x(t),u(t))  x(t0)  = x,  (3.2a) 

In  addition,  the  m-dimensional  input  vector  u(t)  is  bounded  by  limiting  values  Umax  and 
Umin , according  to  the  inequality 

Umax  — U(t)  < Umin  (3.2b) 

where  the  inequality  signs  are  interpreted  elementwise.  The  general  structure  of  the 
control  law  sought  in  the  new  method  proposed  here  is  a suboptimal  time-invariant 
feedback  controller.  A key  motivation  for  this  approach  is  that,  in  contrast  to  optimal 
open-loop  controllers,  time-invariant  feedback  controllers  enhance  robustness  and 
stability  (Goh  and  Edwards,  1994).  This  makes  them  more  attractive  from  the  design 
point  of  view.  In  addition,  optimal  time-varying  feedback  controllers  are  not  always 
convenient  to  deploy  due  to  the  large  amount  of  information  needed  in  implementation 
(Lewis,  1986).  In  contrast,  suboptimal  time-invariant  feedback  control  laws  are  easier  to 
implement  and  therefore  more  attractive  from  a practical  point  of  view. 
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In  this  work  we  directly  design  suboptimal  time-invariant  feedback  controllers  for 
a general  performance  measure  without  first  obtaining  the  optimal  controller.  In  the 
proposed  approach  the  input  u(t)  in  system  (3.1)  and  (3.2)  is  set  as  the  output  of  a 
backpropagation  NN  structure  with  logistic  activation  functions  and  whose  inputs  are  the 
components  of  the  state  vector  x(t).  This  structure  has  the  ability  to  approximate 
arbitrarily  well  any  input -output  mapping  (Homik,  1989).  Furthermore,  it  can  also 
incorporate  the  constraints  on  the  control  variables.  The  weights  of  the  NN  structure  play 
the  role  of  control  parameters  that  must  be  selected  in  a well-defined  optimal  fashion. 
Subsequently,  the  resulting  system  is  expanded  using  Taylor  series  and  then  a higher- 
dimensional linear  approximation  of  the  resulting  polynomic  system  is  obtained  using  a 
minimal-dimension  Carleman  linearization  procedure  (Chapter  2).  The  Carleman 
approximation  of  the  original  system  is  then  integrated  analytically  to  obtain  an  explicit 
expression  for  the  performance  measure  as  a function  of  the  controller  parameters. 
Finally,  the  optimal  controller  parameters  are  obtained  using  standard  optimization 
algorithms  to  yield  the  desired  suboptimal  time-invariant  feedback  control  law. 

A brief  review  of  basic  NNs  concepts  and  the  specification  of  the  suboptimal 
control  law  structure  are  discussed  in  the  next  section. 

3.2  Control  Law  Structure 

The  structure  of  the  classical  backpropagation  NNs  consists  of  a set  of  neurons 
(also  called  nodes  or  cells)  that  are  logically  arranged  into  two  or  more  layers,  namely,  an 
output  and  one  or  more  hidden  layers.  In  addition,  there  is  an  input  layer  that  consists  of 
the  variable  inputs  and  an  extra  input  that  is  set  equal  to  1 (called  bias).  The  inputs  to  the 
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neurons  in  each  layer  come  exclusively  from  outputs  of  the  neurons  in  previous  layers, 
and  outputs  from  these  neurons  pass  exclusively  to  neurons  in  following  layers. 

An  activation  function  f is  applied  to  the  weighted  sums  of  inputs  in  order  to 
produce  each  neuron’s  output.  Thus,  the  neuron’s  output  yj  is  given  by  the  expression 

yrf(ie,,x,+eK.)  (3-3) 

i=l 


where  the  parameters  denoted  0ji  are  the  weights  and  x;  are  the  inputs  to  the  neuron. 

The  networks  usually  have  a single  layer  of  hidden  neurons  in  addition  to  the 
input  and  output  layers.  Such  a configuration  is  called  a three-layer  network.  Rarely 
more  than  one  hidden  layer  is  needed  in  the  network. 

Although  the  operational  characteristics  of  a neuron  are  controlled  by  the  weights, 
the  structure  of  the  activation  function  nevertheless  is  of  relevance.  Basically,  the 
activation  function  is  a nonlinear  operator  that  when  applied  to  the  input  of  a neuron 
determines  the  output  of  that  neuron.  Furthermore,  the  range  of  this  function  is  usually 
between  0 and  1 . The  majority  of  current  network  models  use  a sigmoidal  function.  This 
function  is  continuous,  its  derivative  is  always  positive,  and  its  range  is  bounded.  One  of 
the  most  commonly  used  sigmoidal  forms  is  the  logistic  function 


l + e 

One  major  advantage 
the  form 


of  this  function  is  that  its  derivative 


(3.4) 


can  be  immediately  written  in 


f'(x)  = f(x)(l-f(x))  (3.5) 

In  this  work,  the  control  law  is  designed  with  the  structure  shown  in  Figure  3-1 
that  includes  one  hidden  layer  and  one  output  node  plus  the  logistic  activation  function 
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(3.4).  The  variable  u,  shown  in  the  figure  represents  the  scaled  control  variable  u, . A 
small  number  of  hidden  nodes  should  be  employed  to  keep  the  number  of  parameters 
small.  Often  a single  hidden  node  suffices.  Under  this  design  each  element  of  the  control 
vector  is  given  by  the  control  law 


Ui(x,0)  = 


u max,i  u min,i 


■ + u 


"(ei.(v+ixo + £ e^ixjy  ti ) 

1 + e j=1 


nun.i 


, i = l,2,...,m 


(3.6) 


where 


y*= <3-7> 

-(e^  + Ze^Xj) 

1 + e 

and  n is  the  dimension  of  the  state  vector  and  v is  the  number  of  nodes.  The  vector 
parameter  0 is  of  dimension  m(nv  +2v+l). 

The  structure  (3.6)  incorporates  in  one  analytical  representation  the  nonlinearity 
of  the  relationship  between  the  inputs  and  outputs  as  well  as  the  constraints.  The  logistic 
function  has  a derivative  that  can  be  obtained  easily;  this  facilitates  the  algebraic  steps 
used  later  in  the  linearization  and  optimization  procedure. 

As  a result  of  the  structure  of  Figure  3-1,  the  control  law  is  a function  of  the  states 
x and  of  the  parameters  0,  which  for  brevity  is  denoted  as 

u(t)  = E(x(t),0).  (3.8) 

This  then  changes  system  (3.1)-(3.2)  to 


- M(ia,))+}o(i(t),e)di 


(3.9) 


i(t)  = (p(x(t),0) 


(3.10) 


where 
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Q(i(t),e)=N(i(t),e(i(t),8)) 
cp(i(t),e)  = f(i(t),8(i(t),e)). 


Note  that  the  problem  (3.9)-(3. 10)  implicitly  satisfies  the  constraints  (3.2b)  because  of  the 
form  (3.6)  selected  for  the  control  law. 

3.3  Formulation  of  the  Suboptimal  Problem 

Under  the  assumption  that  the  system  (3.10)  is  analytic  at  the  point  x(to),  a Taylor 
series  expansion  with  respect  to  i gives 


/th  Kronecker  product,  i.e.,  x'[2]-x0x . Note  for  example  that  if  x = [x,  x2]T  then 


states,  namely,  one  of  the  monomials  x,x2  and  x2x,  . 

Let  the  vector  of  monomials  that  includes  the  first  k Kronecker-product  vectors  be 
defined  as 


Differentiating  x'  with  respect  to  time  yields  a higher-dimensional  equivalent  to  (3. 1 1): 


*=a10+2v'1" 


(3.11) 


x'[2)  ={x,2  x,x2  x2x,  x2]t,  hence,  the  second  Kronecker  product  vector  has  redundant 


(3.12) 


(3.13) 


where 
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a;, 

A' 

Ajjt.i 

A ' 

" lk 

a;0 

A' 

"21 

A2,k-2 

A 2,k-l 

K = 

0 

A' 

"30 

A 3,1c  .3 

A 3tk-2 

0 

0 ••• 

Aw 

a;, 

with  A'j  = A(j  ®Ini-, 

+ I„  ®A'_ 

jfori 

= 2,  3, 

, k,  andj  = 0,  1, 


(3.14) 


k-i+1,  and  where  I n is 


the  man  identity  matrix. 

Lumping  into  single  terms  all  components  of  x'  that  are  repetitions  of  the  same 
monomial,  the  redundant  states  are  eliminated.  For  example,  in  the  two-dimensional  case 
x'  = [x,  x2  x,2  x,x2  x2x,  x2]T  is  reduced  to  x'2  = [x,  x2  x2  x,x2  x2]T.  The 
matrix  S'k  must  then  be  modified  by  adding  columns  that  correspond  to  permutations  of 
the  same  monomial  to  yield  one  column,  and  by  omitting  the  corresponding  rows. 
Corresponding  rows  of  g'k  must  also  be  omitted.  This  does  not  modify  the  structure  of 
the  matrices  defined  in  (3.14)  but  reduces  the  size  of  the  submatrices  A-j  and  (except 

for  AJ,  and  A,'0 ) as  well  as  the  size  of  the  zero  matrices.  Hence,  matrix  Sk  and  vector  gk 
are  respectively  replaced  by  Sk  and  gk , of  reduced  dimensions,  obtained  by  replacing 
A'j  in  (3.14)  by  the  reduced  AtJ . 

Truncating  all  terms  of  order  higher  than  Ar  in  (3.13)  yields  the  approximate  linear 

system 


xk  =Sk(0)xk +gk(0)  (3.15) 

where  xk  approximates  the  state  vector  xk  . This  linearization  procedure  was  first 
introduced  by  Carleman  (1932)  and  has  been  used  for  several  nonlinear  problems  (cf. 


Krener,  1974;  Svoronos  et  al.,  1980;  Tsiligiannis  and  Lyberatos,  1987;  Kowalski  and 
Steeb,  1991).  The  integration  of  (3. 1 5)  yields 


(3.16) 


The  Taylor  series  expansion  of  the  first  term  and  of  the  integrand  of  the 
performance  measure  (3.9)  at  x(t)  = x(to)  and  subsequent  truncation  of  all  terms  of  order 
£+1  or  higher  yields 


The  performance  measure  (3.18)  is  an  algebraic  expression  that  can  be  optimized  for  the 
weights  (0).  Therefore,  the  dynamic  optimization  problem  (3.1)-(3.2)  is  transformed  to 
an  unconstrained  algebraic  optimization  problem.  The  suboptimal  feedback  control  law 
u(x,  0° ) is  of  the  form  (3.6)  where  0°  = arg  minJ(0) . 

3.4  Example 

A one-dimensional  example  presented  by  Kirk  (1970)  is  used  to  evaluate  the  new 
approach.  Consider  the  minimization  of  the  performance  measure 


(3.17) 


Substituting  (3.16)  in  (3.17)  and  integrating  yields 


j(0)  = c,  - cT|(I  - es‘<6x"--))S1'(e)gl  (0)]+  bo(0)(tf  - 10) 
-bT(0)[si(0)(tt-to)  + I-es'™"->]sk!(0)g1(e) 


(3.18) 


(3.19) 


0 


subject  to 
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x = ax(t)  + u(t)  , x(0)  = x0  (3.20) 

where  H = 5,  a = 0. 1,  tf=  15,  and  Xo  = 1 . We  assume  that  the  constraints  on  the  input  are 
Umin  ~0.5  and  u^x  0.5. 

The  control  law  for  this  model  is  of  the  form  (3.6),  i.e., 

uRe)=l  + el.e.,»-a5  <32» 

where 


y i + e-(e.o+e„x)  (3.22) 

and  where  the  parameter  vector  is  0 = [01O  0„  02O  02]  J7  . 

Changing  (3.19)  and  (3.20)  to  deviation  variables  and  substituting  the  control  law 

yields 


x = 0.  lx(t)  + 


1 


\ + e-(eM+0jiy) 


-0.4  , x(t0)  = 0 


and 


(3.23) 


J = f(x(15)  + 1)’+I|[T 


1 


12 


+ e 


-(e2o+eJiy) 


-0.5 


dt 


(3.24) 


o 

In  this  example  we  design  three  suboptimal  control  laws  using  Carleman 
linearization  approximations  of  second-  ( k = 2),  fourth-  (k  = 4),  and  sixth-  (k  = 6)  order. 
The  second-order  Carleman  linearization  of  (3.23)  is  of  the  form  (3.15)  with  k = 2 
yielding 


S2  = 


'a, 

V 

2 

5 §2 

0 

_2f0  2a,  _ 

where 


a,  = o.i +uy 

a2  =u”(yY+u'y 
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f0  -0.4+  _ .(e20+e21y|1W)) 


1 + e 


yL=7 


+ e 


■®10 


and  where  the  prime  notation  indicates  the  following  derivatives: 


u 


f 


du(x,  0) 

ay 


x=0 


n •(02o+02iy|!„o) 

d2,c 


a2u(x,e> 

e^e  <e2,+e2‘y|,'°) 

g-c^o+^iylj.o)  _ |1 

ay2 

A second-order  Taylor  series  expansion  of  (3.24)  yields  the  parameters 


5 


1 + e-<62o+92.y|1.0) 


0.5 
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where 

b,  = z'y ' 

b2  = z'(yY  +z'y* 

and  where 


z = 


1 

1 + e-(ett+6iiy) 


z = 


dz 

dy 


x=0 


q e-(ew+e3iy|».o) |j  _ j 

jj  _l_  g-^o+^iylj-o)  J3 


d2z 

q2  _-(e20+62iy|».o) 
_ _y2!e 

jl  _ ^g^^o+^yl*. 

o)  _|_  g’2(62O+02iyk.o)  1 

dy2 

x=0  [ 

_|_  g-(®20+®2iy|«-o) 

The  approximate  performance  measure  J(0),  as  given  by  (3.18),  is  then 
minimized  using  the  simplex  algorithm  ofNelder  and  Mead  (1964)  to  yield 
0°o  =0.41168 

0°  = 0.02250 
0°o  =-0.14818 
0°,  = -0.62974 

and  J(0°),  according  to  equation  (3.18),  was  0.0620  as  shown  in  Table  3-1. 

A fourth-order  Carleman  linearization  of  (3.23)  yields  an  equation  of  the  form 


(3.15)  with  k = 4 
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a, 

a^ 

2 

6 

24 

S4  = 

2f0 

2a, 

a2 

3 

0 

3f0 

3a, 

3a2 

2 

0 

0 

4f0 

4a, 

where 

a,  = 0.1  + u'y' 

a2  =u'(y,)2+u'y' 

a 3 =u*(y')3+3u'y'y'  + u'y' 


a4  = u'"’(y')4  +6uV(y')2  + 4uyy"  + 3u'(y'')2  +u'y,m 


The  performance  measure  (3.24)  expansion  coefficients  are 


C0 


5 

2 


c 


T 


5 

2 


0 0 


b0 


4 


1 

2 _i_  e^e»+e2i  y|f.0) 


bT 


h.  h. 

2 6 24 


where 


b,  = z'y' 

b2=z'(y')2+z'y' 

b3  = z"(y')3 +3z’y  y'  + z'y" 


b4  = z'"(y')4  + 6z*y'(y')2  +4zVy',  + 3z'(y')2  + z'y" 
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The  optimal  parameter  values  are 
e,°0  =-0.19273 

0°  =0.36337 
0°o  = 2.46740 
0°,  = -7.25665 

and  the  corresponding  calculated  J(0)  is  given  in  Table  3-1. 

Details  of  the  sixth-order  system  are  omitted  for  brevity  but  the  results  are  also 
presented  in  Table  3-1. 

The  time- varying  optimal  control  law  for  the  problem  (3.19)  and  (3.20)  is  given 
by  Kirk  (1970)  as 

u(t)  = -2K(t)x(t)  (3.25a) 


A simulation  study  was  carried  out  to  evaluate  the  performance  of  each  of  the 
three  suboptimal  time-invariant  feedback  solutions  with  respect  to  the  performance  of  the 


variable  prescribed  by  each  of  the  four  control  laws  investigated.  Note  that  the  second- 
order  suboptimal  law  fails  to  approximate  the  optimal  solution;  however,  the  higher  order 


converge  towards  the  time-varying  optimal  input  trajectory.  A similar  behavior  is 
observed  on  Figure  3-3,  where  the  higher-order  suboptimal  control  laws  are  shown  to 
produce  state  trajectories  that  converge  towards  the  optimal  trajectory. 


where 


(3.25b) 


optimal  time- varying  solution  (3  .25).  Figure  3-2  shows  the  time  evolution  of  the  control 


suboptimal  control  laws  based  on  fourth  and  sixth-order  approximations  progressively 
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The  level  of  suboptimality  of  the  three  feedback  control  laws  designed  is 
quantified  in  Table  3-1  where  the  performance  index  (3.24)  is  calculated  from  simulation 
results  using  gaussian  quadrature  to  evaluate  the  integral  term.  The  table  shows  that,  as 
the  order  of  the  Carleman  linearization  increases,  the  time-invariant  feedback  control  law 
(3.21)  converges  to  a controller  with  performance  measure  very  close  to  its  optimal 
value.  The  table  also  shows  that  the  approximate  performance  functional  given  by  (3.18) 
gives  values  close  to  those  obtained  by  integration. 
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Table  3-1 . Performance  measure  for  the  example 


Method 

J 

(by  integration) 

J 

(by  eq.(3.l8)) 

2nd  order  Carleman 

0.0620 

0.0620 

4th  order  Carleman 

0.0525 

0.0511 

6th  order  Carleman 

0.0524 

0.0523 

Time-varying  optimal  control 

0.0523 

— 
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ui  ’ umin,i 

ui  = 

umax,i  ~ u min,i 


Figure  3-1.  Structure  of  proposed  control  law  featuring 
a bias  term,  one  hidden  layer  and  one  output  layer. 
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Figure  3-2.  Evolution  of  the  control  variable  in  the 
example  as  a function  of  time  for  the  optimal  control 
law  and  for  three  suboptimal  control  laws  of  varying 
order. 
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Figure  3-3.  Evolution  of  the  state  in  the  example  as 
a function  of  time  using  the  optimal  control  law  and 
three  suboptimal  control  laws  of  varying  order. 


CHAPTER  4 
CONCLUSIONS 


Two  new  methods  had  been  proposed  in  this  work.  The  first  method  which  is 
used  for  discretizing  nonlinear  continuous-time  systems  is  proposed  in  Chapter  2.  The 
method  is  based  on  the  modification  of  the  Carleman  linearization.  The  technique 
involves  an  optimization  of  a performance  measure  which  is  used  to  determine  the 
unknown  coefficients.  The  proposed  discretization  method  yields  models  that  are  more 
accurate  than  both  the  classical  and  other  more  recently  developed  methods.  Moreover, 
the  performance  of  the  low  order  versions  of  the  new  method  is  considerably  accurate 
when  compared  to  the  higher  order  versions  of  the  Carleman  method,  especially  for  large 
sampling  periods.  Furthermore,  the  use  of  different  cases  shows  the  trade  off  between 
simplicity  and  accuracy. 

The  second  method,  proposed  in  Chapter  3,  is  used  for  the  design  of  near-optimal 
feedback  control  laws.  The  new  approach  is  based  on  using  Carleman  linearization 
coupled  with  a concept  similar  to  that  used  in  backpropagation  neural  networks.  The  key 
idea  for  the  synthesis  of  the  control  law  utilizes  the  fact  that  Backpropagation  Neural 
Networks  can  represent  arbitrarily  and  accurately  any  nonlinear  relationship  between 
inputs  and  ouputs  using  a form  of  sigmoidal  functions.  This  function  can  be  explicitly 
stated  by  a set  of  parameters  (weights)  that  are  obtained  by  optimizing  an  algebraic 
performance  measure.  The  proposed  method  yields  a near  optimal  control  law  that  can 
adequately  approximate  the  optimal  control  law. 
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New  avenues  are  proposed  for  future  research.  The  effect  of  increasing  the 
number  of  nodes  used  in  the  network  on  the  accuracy  of  the  results  can  be  studied.  Other 
forms  of  the  sigmoidal  function,  such  as  the  hyperbolic  tangent  or  the  modified 
hyperbolic  tangent,  can  be  used  to  examine  if  more  complex  models  will  lead  to  better 
performance.  This  will  certainly  affect  the  optimization  problem  which  in  turn  may  or 
may  not  justify  the  increasing  model  complexity.  Furthermore,  it  may  be  of  significant 
interest  to  seek  the  development  of  the  conditions  for  convergence,  stability,  and  extent 
of  suboptimality  for  the  new  method.  Another  avenue  is  to  use  a modified  Carleman 
linearization  (along  the  lines  of  Chapter  2)  for  the  design  of  suboptimal  time-invariant 
feedback  controllers. 


APPENDIX  A 

MATLAB  CODES  FOR  EXAMPLE  ONE  SECTION  2.5  (FIRST  ORDER) 


function  [ t, x, u] =closed (T) 

x ( 1 ) = 0; 

t (1 ) = 0; 

alpha  = 1; 
uO  =0.001; 

for  i=l:10 


r = i; 

if (r<=5) 
xsp= . 05*r; 
else 

xsp= . 05* ( 10-r) ; 
end 


xl  = x (i) ; 
par(l)=  xsp; 
par (2 ) = xl; 
par (3) = T; 
par(4)=  alpha; 

u=  f solve (' ddsl uO, [], par) ; 


vl 

v2 

s = 

a 

b 

f 

c 

t (i+l)  = 
x (i+1 ) = 


l+2*alpha+u; 
alpha* u; 

sqrt (vlA2+4*alpha*v2) ; 
exp (s*T) ; 

2 * alpha* x (i) +vl; 

(b+s) *a/ (b-s) ; 

(f+1)/ (f-1) ; 
i; 

(-vl+c*s) / (2*alpha) ; 


end 
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function  f f=ddsl (u, par ) 

n = l+2*par (4) +u; 
m = exp (- (n) *par (3) ) ; 

ff  = par (1) -m*par (2) -par (4) *u* (1-m) /n; 


Open  Loop  Response  for  first  order  Carleman 
T=0 . 05 

Time  x 


1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


0.04088430967938 

-0.00590891515539 

0.03602224949802 

-0.01028634716696 

0.03242034672561 

-0.01352922856093 

0.02975199152269 

-0.01593161418509 

0.02777522536911 

-0.01771134522857 


Time 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4 . 00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


X 

0.19793675271080 

-0.20933256581380 

0.16815415625682 

-0.21975462057959 

0.16667136809100 

-0.22027350413275 

0.16659754441521 

-0.22029933782369 

0.16659386895082 

-0.22030062400742 
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T=5 


Time 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


X 

0.23076922998502 

-0.42855327267078 

0.23076922852868 

-0.42855327267082 

0.23076922852868 

-0.42855327267082 

0.23076922852868 

-0.42855327267082 

0.23076922852868 

-0.42855327267082 


Closed  Loop  Response  for  first  order  Carleman 


T=0 . 05 

Time 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


T=0 . 3 

Time 

1 . 00000000000000 

2 . 00000000000000 

3.00000000000000 

4 . 00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 


X 

0.04995840963645 

0.09972646630620 

0.14927638007948 

0.19861590036399 

0.24775377770864 

0.19769214366122 

0.14860157796435 

0.09927685528898 

0.04973386008177 

-0.00003711777452 


x 

0.04975918863798 

0.09869254823517 

0.14681253312266 

0.19424165617951 

0.24108351479316 

0.19211875797553 

0.14514677748533 

0.09749054864362 

0.04911105401885 
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10.00000000000000 


Time 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


-0.00009434621061 


x 

0.04923217582599 

0.09716753409927 

0.14411518220418 

0.19033709383364 

0.23606796755289 

0.19033709333345 

0.14411518080360 

0.09716753097860 

0.04923217098373 

-0.00000011136524 


global  x J 

J1  = 24; 
gama=  [ - 1 1 ] ; 

optgama  = fmin ( ' discret9 ' , gama ( 1 ) , gama (2 ) ) 

if ( J<J1 ) J1=J 
Min= J1 ; 
end 

function  J=discret9 (gama) 
global  x J 
T = 5; 
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alpha=  1 ; 

u = [-.9  -.8  -.7  -.6  -.5  -.4  -.3  -.2  -.1  le-04  .1  .2  .3 
.4  .5  .6  .7  .8  .9  1 1.1  1.2]; 

[n, m] = size (u) ; 

for  j=l:m 

x(l)  = 0; 

for  i=l : 1 


vl 

= l+2*alpha+u ( j ) ; 

v2 

= alpha*u ( j ) ; 

gam 

= gama ( 1 ) *u ( j ) +gama (2 ) *u ( j ) *2 ; 

LI 

= - ( vl+alpha*gam) ; 

a 

= exp (L1*T) ; 

x (i+1) 

= a*x (i) +v2* ( 1-a) / (-L1 ) ; 

end 

w(j)  = ( (x  (i+1)  A2-gam*x  (i+1)  ) /x  (i+1)  /'2)  A2; 
end 

J = 0; 

for  i=l:m 

J = J+w (i) ; 

End 

J 


function  [ t, x] =discreta (T) 


t (1) 
x(l) 

alpha  = 
gama ( 1 ) = 
gama ( 2 ) = 
gama (3) = 


0; 

0; 

l; 

0.04643071097378; 

-0.00118452449519; 

0.00002158383575; 


for  i=l : 10 
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r 

ur 

vl 

v2 

% gam 
% gam 
gam 
LI 
A 


i-l; 

0.9* (-1) Ar; 
l+2*alpha+ur; 
alpha*ur; 
gama (1) *ur; 

gama ( 1 ) *ur+gama (2 ) *urA2 ; 
gama (1) *ur+gama (2) *urA2+gama (3) *urA3; 
- (vl+alpha*gam) ; 
exp (L1*T) ; 


t(i+l)  = i; 

x (i+1)  = a*x (i) +v2* (1-a) / (-L1) ; 


end 


function  [ t, x] =closedb (T) 

x (1)  = 0; 

t (1 ) = 0; 

alpha  = 1; 
uO  = 0.001; 

for  i=l:10 


r = i; 


if (r<=5) 
xsp= . 05*r; 
else 

xsp=. 05* (10-r) ; 
end 

xl  = x (i) ; 

par(l)=  xsp; 
par (2 ) = xl; 
par (3) = T; 
par(4)=  alpha; 

u = f solve ( 'dds3' , uO, [] ,par) ; 

vl  = l+2*alpha+u; 

v2  = alpha*u; 

s = sqrt (vl A2+4*alpha*v2 ) ; 
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a = exp (s*T) ; 

b = 2*alpha*x (i) +vl; 

f = (b+s) *a/ (b-s) ; 

c = (f+1) / ( f — 1 ) ; 

t (i+1 ) = i; 

x (i+1 ) = (-vl+c*s) / (2* alpha) ; 


end 


function  f f=dds3 (u, par) 

gama(l)  = 0.04643071097378; 

gama ( 2 ) = -0.00118452449519; 

gama ( 3 ) = 0.00002158383575; 

vl  = l+2*par (4) +u; 

v2  = par (4 ) *u; 

%gam  = gama(l)*u; 

%gam  = gama  ( 1 ) *u+gama  (2 ) *u/'2 ; 

gam  = gama (1) *u+gama (2) *uA2+gama (3) *uA3; 

LI  =- (vl+par (4 ) *gam) ; 

A = exp (Ll*par (3) ) ; 

ff  = par (1) -a*par (2) -v2* (1-a) / (-L1) ; 


First  order  linear  in  u 
T=0 . 05 

optgama  = 0.04623651268079 
Min  = 0.00575293637520 


Open-loop  response 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4 . 00000000000000 

5.00000000000000 


0.04084318634951 

-0.00591304317091 

0.03598784224882 

-0.01029353329161 

0.03239091487136 
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6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


Closed-loop 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


First  order  for  u/'2 
T=0 . 05 


-0.01353868018853 

0.02972624553166 

-0.01594274413856 

0.02775220993273 

-0.01772371851643 


0.05002014959766 

0.09993833280848 

0.14968910607864 

0.19928951441184 

0.24875963763025 

0.19749044793049 

0.14840775659114 

0.09909449409661 

0.04960210921308 

-0.00008651223837 


optgama  = 0.04643885281353 
-0.00117768214738 

Min  = 6. 171336653536668 e- 007 


Open-loop  response 

1.00000000000000 

2 . 00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8 . 00000000000000 

9.00000000000000 

10.00000000000000 


0.04084394857585 

-0.00591145599096 

0.03598972042235 

-0.01029118807240 

0.03239327670237 

-0.01353608292312 

0.02972871149820 

-0.01594018947492 

0.02775456473543 

-0.01772136531524 


Closed-loop 
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1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


0.05001867584492 

0.09993204421170 

0.14967414088016 

0.19926001674926 

0.24870686944350 

0.19748888696379 

0.14840420296765 

0.09909050303910 

0.04959877087704 

-0.00008791055794 


First  order  for  uA3 
T=0 . 05 

optgama  = 0.04643071097378 
-0.00118452449519 
0.00002158383575 

Min  = 4. 47447295963835 9e-011 


Open-loop  response 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8 . 00000000000000 

9.00000000000000 

10.00000000000000 

Closed-loop 

1 . 00000000000000 

2.00000000000000 

3.00000000000000 

4 . 00000000000000 

5.00000000000000 

6.00000000000000 

7 . 00000000000000 

8 . 00000000000000 

9.00000000000000 


0.04084394574675 

-0.00591144749813 

0.03598972526247 

-0.01029117571471 

0.03239328523128 

-0.01353606950040 

0.02972872128331 

-0.01594017659558 

0.02775457435711 

-0.01772135382565 


0.05001869022451 

0.09993214120996 

0.14967445883637 

0.19926082724944 

0.24870868014239 

0.19748885131141 

0.14840419611289 

0.09909047634033 

0.04959873621645 
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10.00000000000000  -0.00008792923287 


First  order  ACS1.1  (for  linear  u) 
see=gamal*u+ . . . 


T=5 


optgama  = 0.29010839374162 
Min  = 1.35924679647996 


Open-loop  response 


Time 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


Closed-loop  response 
Time 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4 . 00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8 . 00000000000000 

9.00000000000000 

10.00000000000000 


X 

0.21628908897519 

-0.48935067668689 

0.21628908852447 

-0.48935067668693 

0.21628908852447 

-0.48935067668693 

0.21628908852447 

-0.48935067668693 

0.21628908852447 

-0.48935067668693 


x 

0.04995202963372 

0.10001422836719 

0.15047490612127 

0.20164564185820 

0.25385835437216 

0.20164564176869 

0.15047490586355 

0.10001422711848 

0.04995202707556 

-0.00000000024856 
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First  order  for  uA2 
T=5 

optgama  = 0.35502884367777 
-0.12862934915581 

Min  = 0.14661197864843 


Open-loop  response 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 

Closed-loop 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7 . 00000000000000 

8 . 00000000000000 

9.00000000000000 

10.00000000000000 


0.21869416207396 

-0.53672868067917 

0.21869416145250 

-0.53672868067931 

0.21869416145250 

-0.53672868067931 

0.21869416145250 

-0.53672868067931 

0.21869416145250 

-0.53672868067931 


0.05005646925102 

0.10020419163971 

0.15030644190391 

0.20007844450633 

0.24908355245331 

0.20007844441330 

0.15030644147426 

0.10020419044178 

0.05005646706795 

-0.00000015079648 


First  order  ACS1.1  (for  linear  u) 
see=gamal*u+ . . . 


T=5 


optgama  = 0.29010839374162 
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Min  = 1.35924679647996 


Open-loop  response 
Time 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


Closed-loop  response 
Time 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


X 

0.21628908897519 

-0.48935067668689 

0.21628908852447 

-0.48935067668693 

0.21628908852447 

-0.48935067668693 

0.21628908852447 

-0.48935067668693 

0.21628908852447 

-0.48935067668693 


x 

0.04995202963372 

0.10001422836719 

0.15047490612127 

0.20164564185820 

0.25385835437216 

0.20164564176869 

0.15047490586355 

0.10001422711848 

0.04995202707556 

-0.00000000024856 


APPENDIX  B 

MATLAB  CODES  FOR  EXAMPLE  ONE  SECTION  2.5  (SECOND  ORDER) 


function  [ t, x, ul ] =closeda (T) 

x(l)  = 0; 
t(l)  = 0; 
ul(l)  = 0; 
alpha  = 1 ; 
uO  = 0.001; 

for  i=l : 10 


r = i; 


if  (r<=5) 
xsp=. 05*r; 
else 

xsp=. 05* (10-r) ; 
end 


xl  = x (i)  ; 

par(l)=  xsp; 
par(2)=  xl; 
par (3) = T; 
par(4)=  alpha; 

u = f solve (' dds2 uO,  [],  par) ; 

vl  = l+2*alpha+u; 

v2  = alpha*u; 

s = sqrt (vlA2+4*alpha*v2)  ; 

a = exp (s*T) ; 

b = 2*alpha*x (i) +vl; 

f = (b+s) *a/ (b-s) ; 

c =(f+l) / (f-1) ; 

t (i+1)  = i; 

x(i+l)  = (-vl+c*s) / (2*alpha) ; 
ul (i+1) = u; 
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end 

function  f f=dds2 (u, par) 

vl  = l+2*par (4) +u; 
v2  = par (4) *u; 

s = sqrt (vlA2-8*par (4) *v2)  ; 

LI  = (-3*vl+s) /2; 

L2  =(-3*vl-s) /2; 
a = exp (Ll*par (3) ) ; 
b = exp (L2*par (3) ) ; 

P = [a  0;  0*b]  ; 

E = [par (4)  par (4) ;-vl-Ll  -vl-L2]; 

Einv=  inv(E); 

R = [ (a-1) /LI  0; 0 (b-1) /L2] ; 

C = E*P*Einv; 

D = E*R*Einv; 

ff  = par  (1)  -c  (1, 1)  *par  (2)  -c  (1, 2)  *par  (2)  /'2-v2*d  (1,1)  ; 


Open  Loop  Response  for  second  order  Carleman 
T=0 . 05 


Time 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


X 

0.04085649221859 

-0.00595578406852 

0.03596488126783 

-0.01035332567148 

0.03235138713174 

-0.01360315892526 

0.02968004560156 

-0.01600637549477 

0.02770410997889 

-0.01778437737920 


T=0 . 5 


Time  x 

1.00000000000000  0.19222315004222 

2 . 00000000000000  -0.21698158629471 
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3.00000000000000 

4 . 00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


Time 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


0.16384476570641 

-0.22741338931106 

0.16241131983133 

-0.22794529717256 

0.16233805990464 

-0.22797249454516 

0.16233431355445 

-0.22797388539209 


x 

0.21787709491554 

-0.53798049421949 

0.21787709474172 

-0.53798049421958 

0.21787709474172 

-0.53798049421958 

0.21787709474172 

-0.53798049421958 

0.21787709474172 

-0.53798049421958 


closed  Loop  Response  for  second  order  Carleman 
T=0 . 5 


Time 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


X 

0.05001150364399 

0.10008689209373 

0.15025584667553 

0.20052827088905 

0.25088676400700 

0.20052827091270 

0.15025584676422 

0.10008689227392 

0.05001150384656 

-0.00000014813352 
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T=5 


Time 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


global  x J 

J1  = 100; 

gama  = [ - 1 1 ] ; 


x 

0.05001150364399 

0.10008689209373 

0.15025584667553 

0.20052827088905 

0.25088676400700 

0.20052827091270 

0.15025584676422 

0.10008689227392 

0.05001150384656 

-0.00000014813352 


optgama  = fmin ( ' discnw6 ' , gama ( 1 ) , gama (2 ) ) 

if ( J< J1 ) J1=J; 

Min=Jl ; 

End 


Min 


function  J=discnw6 (gama) 

global  x J 

T = 0.05; 
alpha  = 1; 

u =[-.9  -.8  -.7  -.6  -.5  -.4  -.3  -.2  -.1  le-04  .1  .2  .3 
.4  .5  .6  .7  .8  .9  1 1.1  1.2]; 

[n, m]  = size (u) ; 

for  j=l:m 

x (1)  = 0; 
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for  i=l : 1 


% 

% 

% 

% 


% 


vl  = l+2*alpha+u ( j ) ; 

v2  = alpha*u(j); 

seel=  gama (1) *u ( j ) ; 

seel=  gama (1) *u ( j ) +gama (2) *u ( j ) A2; 

seel=  gama (1) *u ( j ) +gama (2) *u ( j ) A2+gama (3) *u ( j ) A3; 

seel=  gama (1) *u ( j ) +gama (2) *u ( j ) A2 

+gama (3) *u(j) A3+gama (4 ) *u ( j ) A4 ; 
seel=  gama (1) *u ( j ) +gama (2) *u (j ) A2 

+gama (3) *u ( j ) A3+gama (4 ) *u  ( j ) A4+gama (5) *u ( j ) A5; 
seel=  gama (1) *u ( j ) +gama (2) *u ( j ) A2 

+gama (3) *u ( j ) A3+gama (4) *u ( j ) A4 
+gama (5) *u ( j ) A5+gama (6) *u ( j ) A6; 
seel=  gama (1) *u ( j ) +gama (2) *u ( j ) A2 

+gama (3) *u ( j ) A3+gama (4 ) *u ( j ) A 4+ gama (5) *u  ( j ) A5 
+gama (6) *u ( j ) A6+gama (7) *u ( j ) A7; 


gam(l)=  seelA2; 
gam ( 2 ) = 0 ; 

s = sqrt (4*alphaA2*gam  (2) A2+4*vl*alpha*gam (2 ) 
+8*alphaA2*gam ( 1 ) +vl A2-8*alpha*v2 ) ; 


LI 

= • 

-alpha*gam(2) + (-3*vl+s) / 2; 

L2 

= * 

-alpha*gam(2) + (-3*vl-s) /2; 

A 

= 

exp  (L1*T) ; 

B 

= 

exp (L2*T) ; 

P 

= 

[a  0;  0 b]  ; 

E 

= 

[alpha  alpha; -vl-Ll  -vl-L2] 

Einv 

= 

inv (E) ; 

R 

= 

[ (a-l)/Ll  0; 0 (b-l)/L2] ; 

C 

= 

E*P*Einv; 

d 

= 

E*R*Einv; 

t ( i ) = i ; 

x (i+1)  = c(l, 1) *x(i)+c(l,2) *x (i) A2+v2*d (1,1); 
end 

w ( j ) = ( (x (i+1) A3-gam(l) *x (i+1) 

-gam (2) *x (i+1 ) A2 ) /x ( i+1 ) A3) A2 ; 

end 


J=0; 


for  i=l:m 
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J=J+w (i) ; 
End 
J 

%gama ( 1 ) ; 
%gama (2) ; 


function  [ t, x] =discre2 (T) 


t(l)  = 0; 
x(l)  = 0; 
alpha  = 1; 

gama ( 1 ) = 0.04622966868911; 


for  i=l:10 


r 

ur 

vl 

v2 

seel 
% seel 
% seel 
% seel 

% seel 

% seel 

% seel 


gam ( 1 ) 
gam ( 2 ) 
s 


= i-l; 

= 0.9* (-1) Ar; 

= l+2*alpha+ur; 

= alpha*ur; 

= gama(l)*ur; 

= gama (1) *ur+gama (2) *urA2; 

= gama (1) *ur+gama (2) *urA2+gama (3) *urA3; 

= gama (1) *ur+gama (2) *urA2+gama (3) *urA3 
+gama (4) *urA4; 

= gama (1) *ur+gama (2) *urA2+gama (3) *urA3 
+gama (4) *urA4+gama (5) *urA5; 

= gama (1) *ur+gama (2) *urA2+gama (3) *urA3 
+gama (4) *urA4+gama (5) *urA5 
+gama (6) *urA6; 

= gama (1) *ur+gama (2) *urA2+gama (3) *urA3 
+gama (4 ) *urA4+gama (5) *urA5+gama (6) *urA6 
+gama (7) *urA7; 

= seelA2; 

= 0; 

= sqrt (4*alphaA2*gam (2) A2+4*vl* alpha* gam (2 

4-P*  a 1 r-iH  = A O * rr^Tn  I 1 \ -LttI  \ . 


LI 

=-alpha*gam ( 

L2 

=-alpha*gam ( 

a 

= exp (L1*T) ; 

b 

= exp(L2*T); 

P 

= [a  0;  0 b]  ; 

= [alpha  alpha;-vl-Ll  -vl-L2]; 
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Einv  = inv(E); 

R = [ (a-1 ) /LI  0 ; 0 (b-l)/L2]; 
c = E*P*Einv; 

d =E*R*Einv; 

t (i+1 ) = i; 

x (i+1 ) = c(l/l)*x(i)+c(l,2)*x(i)  A2+v2*d  ( 1/ 1 ) ; 
end 


function  f f=dds4 (u, par) 


% 

% 

% 

% 

% 


gama ( 1 ) = 

vl 

v2 

see  = 

see  = 

see  = 

see  = 

see 

see  = 

gam 

s 

LI 

L2 

a = 

b 


0.04622966868911; 
l+2*par (4) +u; 
par (4) *u; 
gama (1 ) *u; 

gama  ( 1 ) *u+gama  (2 ) *uA2 ; 

gama (1) *u+gama (2) *uA2+gama (3) *uA3; 

gama (1) *u+gama (2) *uA2+gama (3) *uA3+gama (4) *uA4; 

gama (1) *u+gama (2) *uA2+gama (3) *uA3 

+gama (4) *uA4+gama (5) *uA5; 

gama (1) *u+gama (2) *uA2+gama (3) *uA3 

+gama (4) *uA4+gama (5) *uA5+gama (6) *uA6; 

seeA2; 

sqrt (vlA2-8*par (4) *v2+8*par (4) A2*gam) ; 
(-3*vl+s)/2; 

(-3*vl-s)/2; 
exp  (Ll*par  (3) ) ; 
exp (L2*par (3) ) ; 


P = [a  0;  0 b]  ; 

E = [par (4)  par (4) ;-vl-Ll  -vl-L2]; 

Einv  = inv(E); 

R = [(a-1) /LI  0; 0 (b-l)/L2]; 

c = E*P*Einv; 

d = E*R*Einv; 


ff  = par (1) -c (1, 1) *par (2) -c (1, 2) *par (2) A2-v2*d (1,1); 


function  [ t, x] =closedc (T) 

x(l)  = 0; 
t (1)  = 0; 
alpha  = 1; 
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uO  = O.OOl; 
for  i=l : 10 
r = i; 


if (r<=5) 
xsp= . 05*r; 
else 

xsp= . 05* ( 10-r) ; 
end 


xl  = x (i) ; 
par(l)=  xsp; 
par (2 ) = xl; 
par (3) = T; 
par(4)=  alpha; 

u=  f solve (' dds4 uO, [], par) ; 

vl  = l+2*alpha+u; 
v2  = alpha*u; 

s = sqrt (vl A2+4*alpha*v2 ) ; 
a = exp (s*T) ; 
b = 2*alpha*x (i) +vl ; 
f = (b+s) *a/  (b-s) ; 
c = (f+1) / (f-1) ; 

t (i+1)  = i; 

x(i+l)  = (-vl+c*s) / (2*alpha) ; 
end 


Second  order  for  u 
T=0 . 05 

optgama  = 0.04622966868911 

Min  = 0.02236807261794 
Open-loop  response 

1.00000000000000  0.04085654571286 

2.00000000000000  -0.00595564071556 
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3.00000000000000 

4 . 00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


Closed-loop 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


0.03596503265107 

-0.01035311267382 

0.03235158111858 

-0.01360292115257 

0.02968024915604 

-0.01600613912030 

0.02770430440057 

-0.01778415662417 


0.04999990505113 

0.10000012231725 

0.15000201179347 

0.20000676552763 

0.25001535035766 

0.20002853273211 

0.15001713317201 

0.10000499214701 

0.05000106781566 

-0.00000001918236 


Second  order  for  uA2 
T=0 . 05 

optgama  = 0.04643838375997 
-0.00116141640383 

Min  = 2.193201990143124e-006 


Open-loop  response 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


0.04085654379434 

-0.00595563720121 

0.03596503433280 

-0.01035310693464 

0.03235158515656 

-0.01360291404087 

0.02968025471299 

-0.01600613119366 

0.02770431091676 

-0.01778414824219 
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Closed-loop 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


0.04999990958924 

0.10000015248436 

0.15000210166808 

0.20000697946736 

0.25001580524556 

0.20002852913188 

0.15001712364887 

0.10000497843865 

0.05000105367964 

-0.00000002729764 


Second  order  for  uA3 
T=0 . 05 

optgama  = 0.04643070571988 
-0.00116784828476 
0.00002034952925 

Min  = 1.3150172351849 92 e-010 


Open-loop  response 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 

Closed-loop 

1 . 00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 


0.04085654380119 

-0.00595563713325 

0.03596503439293 

-0.01035310683071 

0.03235158524441 

-0.01360291392081 

0.02968025481273 

-0.01600613106940 

0.02770431101894 

-0.01778414812053 


0.04999990954701 

0.10000015204122 

0.15000209985105 

0.20000697388663 

0.25001579046973 

0.20002852913208 

0.15001712361062 
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8.00000000000000 

9.00000000000000 

10.00000000000000 


Second  order  for  u 


0.10000497834719 

0.05000105353608 

-0.00000002740281 


T=5 


optgama  = -0.27545020794361 
Min  = 4.30095719499269 


Open-loop  response 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 

Closed-loop 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


0.21871144283363 

-0.54747543096109 

0.21871144261814 

-0.54747543096121 

0.21871144261814 

-0.54747543096121 

0.21871144261814 

-0.54747543096121 

0.21871144261814 

-0.54747543096121 


0.05000188054247 

0.10000862217782 

0.14998839248929 

0.19988769510133 

0.24962710491168 

0.19988769511689 

0.14998839255227 

0.10000862230858 

0.05000188069910 

-0.00000014812790 


Second  order  for  uA2 


T=5 
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optgama  = 0.35425608332670 
-0.13010761411484 


Min  = 0.61244005239382 


Open-loop  response 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 

Closed-loop 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


0.21849498469303 

-0.56672853897686 

0.21849498447979 

-0.56672853897703 

0.21849498447979 

-0.56672853897703 

0.21849498447979 

-0.56672853897703 

0.21849498447979 

-0.56672853897703 


0.04999740881055 

0.09998814934589 

0.14997564617716 

0.19999927729383 

0.25012899409181 

0.19999927731078 

0.14997564623884 

0.09998814946301 

0.04999740894548 

-0.00000014812423 


Second  order  for  u^3 


T=5 


optgama  = 0.33973593439105 
-0.18360368619533 
0.07915598539186 


Min 


0.11354332858997 
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Open-loop  response 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 

Closed-loop 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


0.21850258375146 

-0.58089540735761 

0.21850258353172 

-0.58089540735784 

0.21850258353172 

-0.58089540735784 

0.21850258353172 

-0.58089540735784 

0.21850258353172 

-0.58089540735784 


0.04999913283630 

0.10000258962767 

0.15001353673411 

0.20002695074230 

0.24997718411539 

0.20002695075968 

0.15001353679943 

0.10000258975408 

0.04999913297945 

-0.00000014812499 


Second  order  for  uM 


T=5 


optgama  = 0.32938112123995 
-0.15916087703063 
0.12722270417632 
-0.05907366747052 

Min  = 0.02504632147283 


Open-loop  response 

1 . 00000000000000 

2.00000000000000 

3.00000000000000 

4 . 00000000000000 

5.00000000000000 

6.00000000000000 


0.21854271550000 

-0.58942448515712 

0.21854271527448 

-0.58942448515737 

0.21854271527448 

-0.58942448515737 
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7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 

Closed-loop 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


0.21854271527448 

-0.58942448515737 

0.21854271527448 

-0.58942448515737 


0.04999956001303 

0.10000179093607 

0.14999868251678 

0.19998631317798 

0.24999692551245 

0.19998631319482 

0.14999868258079 

0.10000179106211 

0.04999956015830 

-0.00000014812549 


Second  order  for  uA5 


T=5 

optgama  = 0.33141108542995 
-0.13749791924959 
0.09684394899844 
-0.10605411696941 
0.04964418719068 

Min  = 0.00618070464791 


Open-loop  response 

1 . 00000000000000 

2.00000000000000 

3.00000000000000 

4 . 00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8 . 00000000000000 

9.00000000000000 

10.00000000000000 


0.21851715739094 

-0.59403963863261 

0.21851715716473 

-0.59403963863289 

0.21851715716473 

-0.59403963863289 

0.21851715716473 

-0.59403963863289 

0.21851715716473 

-0.59403963863289 
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Closed-loop 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 


0.04999920242677 

0.09999900770897 

0.14999609638479 

0.20000240307684 

0.25000845630290 

0.20000240309387 

0.14999609644854 

0.09999900783320 

0.04999920257037 

-0.00000014812540 


Second  order  for  uA6 


T=5 

optgama  = 0.33409308179247 
-0.14232105636675 
0.06242651151102 
-0.07075432691257 
0.09849031048462 
-0.04505941169768 

Min  = 0.00160562970688 


Open-loop  response 

1.00000000000000 

2.00000000000000 

3.00000000000000 

4.00000000000000 

5.00000000000000 

6.00000000000000 

7.00000000000000 

8.00000000000000 

9.00000000000000 

10.00000000000000 

Closed-loop 

1 . 00000000000000 

2 . 00000000000000 

3.00000000000000 

4.00000000000000 


0.21852589709572 

-0.59638151632177 

0.21852589686802 

-0.59638151632205 

0.21852589686802 

-0.59638151632205 

0.21852589686802 

-0.59638151632205 

0.21852589686802 

-0.59638151632205 


0.04999910777536 

0.09999969830725 

0.15000046655569 

0.20000193125578 
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5.00000000000000 

6.00000000000000 

7 . 00000000000000 

8 . 00000000000000 

9.00000000000000 

10.00000000000000 


0.24999324534287 

0.20000193127282 

0.15000046661979 

0.09999969843188 

0.04999910791845 

-0.00000014812526 


APPENDIX  C 

MATLAB  CODES  FOR  EXAMPLE  TWO  SECTION  2.5 
(FIRST  AND  SECOND  ORDER) 


function  [t,xl,x2]=discret5 (T) 
cain=  10; 

cas  = 4.286373240038316; 
tin  = 3e02; 

ts  = 4. 0008647429034 92 e+ 002; 
v = le-02; 
f = 0 . 33e-04 ; 
z = 3e05; 

e = 6e04; 

R = 8.314; 

h = le04; 

rho  = le03; 

C = 4; 

xll  = (cain-cas) /cas; 
x2I  = (tin-ts)/ts; 
a = v*z/f; 
b = e/ (R*ts)  ; 
c = h*cas/ (rho*C*ts) ; 

qs  =-rho*C* (h*v*z*exp (-b) *cas/ (rho*C) +f* (tin-ts) ) ; 
d = qs/ (f*rho*C*ts) ; 

A = a*exp(-b); 

xl(l)  = 0; 
x2 (1)  = 0; 
t (1)  = 0; 

for  i=l : 10 

j = i-1; 
u = 0.9*  (-1) Aj; 

Fl=  [ - ( 1+A)  -b*A; c*A  -1+A*b*c]; 

% eig (FI) 

F = [-(1+A)  -b*A  0 -b*A  -A* (bA2-2*b) /2 


80 


81 


c*A  -1+A*b*c 

0 

b*c*A 

A*c*  (b/N2-2*b)  / 2 

0 0 

-2* (1+A) 

-2*b*A 

0 

d*u  0 

C*A 

-2-A+b*c*A 

-b*A 

0 2*d*u 

0 

2*c*A 

2* (-l+b*c*A)  ] ; 

[E,D] 

= eig (F) ; 

% eig(F) 

all  = 

exp (D (1, 1) *T) ; 

a22  = 

exp (D (2, 2 ) *T) ; 

a33  = 

exp  ( D ( 3 , 3 ) * T ) ; 

a44  = 

exp (D (4, 4) *T) ; 

a55  = 

exp (D (5, 5) *T) ; 

Pi  = [all  a22  a33  a44  a55] ; 

PI  = diag (Pi, 0) ; 

Einv=  inv(E); 

Ri  = [ (all-1) /D(l, 1)  (a22-l ) /D (2, 2)  (a33-l) /D (3, 3) 

(a44-  1 ) /D ( 4 , 4 ) (a55-l) /D (5, 5) ] ; 

Rl  = diag (Ri, 0) ; 

cl  = E*Pl*Einv; 

dl  = E*Rl*Einv; 

t (i+1)  = i; 

xl (i+l)=  cl (1, 1) *xl (i) +cl (1,2) *x2 (i)+cl (1, 3) *xl (i) A2 
+cl (1, 4) *xl (i) *x2 (i)+cl (1,5) *x2 (i) A2 
+d*u*dl (1,2) ; 

x2 (i+1) = cl (2, 1) *xl (i) +cl (2,2) *x2 (i) +cl (2, 3) *xl (i) A2 
+cl (2, 4) *xl (i) *x2 (i)+cl (2,5) *x2 (i) A2 
+d*u*dl (2,2) ; 

end 


function  [t, xl, x2] =carl3 (T) 


cain=  10; 

cas  = 4.286373240038316; 
tin  = 3e02; 

ts  = 4. 0008647429034 92 e+ 002; 
v = le-02 ; 
f = 0 . 33e-04 ; 
z = 3e05; 


e = 6e04; 
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R = 8.314; 

h = le04 ; 

rho  = le03; 

C = 4; 


xll  = (cain-cas) /cas; 
x2I  = (tin-ts) /ts; 
a = v*z/f; 
b = e/ (R*ts) ; 
c = h*cas/ (rho*C*ts) ; 

qs  =-rho*C* (h*v*z*exp (-b) *cas/ (rho*C) + f* (tin-ts) ) ; 
d = qs/ (f*rho*C*ts) ; 

A = a*exp(-b); 

xl(l)=  0; 
x2 ( 1 ) = 0; 
t (1)  = 0; 

for  i=l : 10 


j = i-1; 
u = 0.9* (-1) Aj; 


F= [ - ( 1+A)  -b*A  0 -b*A  (b-.5*bA2)*A 
(b-.5*bA2)*A  (6*bA2-bA3-6*b) *A/6 


0 0 


C*A 

-1+A*b*c  0 

b*c*A  A*c* (bA2-2*b) /2 

0 0 

A*c* (bA2-2*b) /2 

c* (bA3-6*bA2+6*b) * 

A/ 6 

0 

0 -2* (1+A) 

-2*b*A  0 0 

-2*b*A 

2*b*A-bA2*A 

0 

d*u 

0 c*A 

-2-A+b*c*A  -b*A 

0 

b*c*A 

(b*c+ 

. 5*c*bA2-b) *A 

(- . 5*bA2+b) *A 

0 

2*d*u  0 2*c*A  2*(-l+b*c*A) 

0 

0 

2*c*b*A 

(b-2) *c*b*A 

0 

0 0 0 

0 -3* (1+A) 

-3*b*A 

0 

0 

0 0 
-2*b*A 


d*u  0 0 

c*A 


0 -3-2*A+b*c*A 


000  2*d*u  0 

-3-A+2*b*c*A  -b*A 


0 


2*c*A 
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0 

] ; 

[E, D]  = eig(F)  ; 

% eig(F) 

all  = exp (D  (1, 1)  *T) ; 
a22  = exp (D (2, 2) *T) ; 
a33  = exp (D (3, 3) *T) ; 
a44  = exp (D (4, 4) *T) ; 
a55  = exp (D (5, 5) *T) ; 
a66  = exp (D ( 6, 6) *T) ; 
all  = exp (D (7, 7) *T) ; 
a88  = exp (D (8, 8) *T) ; 
a99  = exp (D (9, 9) *T) ; 

Pi  = [all  a22  a33  a44  a55  a66  all  a88  a99]; 

PI  = diag (Pi, 0) ; 

Einv=inv (E) ; 

Ri  =[ (all-1) /D(l, 1)  (a22-l ) /D (2, 2)  (a33-l) /D (3, 3) 

(a44-l)/D(4,4)  (a55-l) /D (5, 5)  (a66-l ) /D ( 6, 6) 

(all-1) /D (7, 7)  (a88-l) /D (8, 8)  (a99— 1) /D(9, 9) ] ; 

Rl  = diag (Ri, 0) ; 
cl  = E*Pl*Einv; 
dl  = E*Rl*Einv; 


000  0 3*d*u  0 

3*c*A  -3+3*b*c*A 


t (i+l)=i; 

xl (i+1) =cl (1, 1) *xl (i)+cl (1,2) *x2 (i)+cl (1, 3) *xl (i) A2+cl  (1,4) 
*xl (i) *x2 (i)+cl (1,5) *x2 (i) A2+cl (1, 6) *xl (i) A3+cl (1, 7) *xl (i) A 
2*x2 ( i ) +cl (1, 8) *xl (i) *x2 (i) A2+cl (1,9) *x2 (i) A3+d*u*dl (1,2) ; 
x2 (i+1 ) =cl (2, 1) *xl (i)+cl (2,2) *x2 (i)+cl (2,3) *xl (i) A2+cl (2,4) 
*xl (i) *x2 (i)+cl (2,5) *x2 (i) A2+cl (2, 6) *xl (i) A3+cl (2, 7) *xl (i) A 
2*x2 (i)+cl (2, 8) *xl (i) *x2 (i) A2+cl (2, 9) *x2 (i) A3+d*u*dl (2,2); 
end 


global  x J 
Jl=1000; 

gama= [ -0.11792836252674 

0.05160393362069] ; 

options ( 1 ) =0 ; 
options (2) =le-10; 
options (3) =le-10; 
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options (14) =20000; 

optgama=fmins ( ' discn6b ' , gama, options) 

if ( J<J1)  J1=J; 

Min=Jl; 

end 


function  J=discn6b (gama) 
global  x J 
T= . 5 ; 
cain=10; 

cas=4. 28637324003831 6; 
tin=3e02; 

ts=4 .0008 647429034 92e+002; 

v=le-02; 

f= . 33e-04 ; 

z=3e05; 

e=6e04 ; 

R=  8.314; 
h=le04 ; 
rho=le03; 

C=4 ; 

xli= (cain-cas) /cas; 
x2i= (tin-ts) /ts; 
a=v*z/f ; 
b=e/ (R*ts) ; 
c=h*cas/ (rho*C*ts) ; 

qs=-rho*C* (h*v*z*exp (-b) *cas/ (rho*C) +f* (tin-ts) ) ; 
d=qs/ (f*rho*C*ts) ; 

A=a*exp (-b) ; 

%u= [0.9  -0.9] ; 

u= [ 0 . 9 .8  .7  .6  .5  .4  .3  .2  .1  .0001  -.1  -.2  -.3  -.4  -.5 
.6  -.7  -.8  -0.9] ; 

[n,m] =size (u) ; 
for  j=l:m 
xl (1 ) =0; 
x2 ( 1 ) =0; 
x3  ( 1 ) =0 ; 
x4 (1) =0; 
x5 (1 ) =0; 
for  i=l:l 

seel=gama (1) *u ( j ) ; 
see2=gama (2) *u ( j ) ; 
gam(l)=seel/v2; 
gam (2) =0; 
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gam  ( 11 ) =see2A2 ; 

gam ( 12 ) =0; 

gam (3) =seel*see2/2; 

gam(4)=seelA2/2; 

gam ( 5) =0; 

gam ( 6) =0; 

gam  (7)  =see2A2/2; 

gam (8) =seel*see2/2; 

gam (9) =0; 

gam(10)=0; 

wll=-l-A+  (b-.5*bA2)  *A*gam(7)  ; ; 

wl2=-b*A+ (b- . 5*bA2 ) *A*gam ( 8 ) + (bA2- . 167*bA3-b) *A*gam(ll) ; 
wl3=0; 

wl4=A* (-b+ (b-.5*bA2) *gam(10) ) ; 

wl5=-  ( . 5*bA2-b)  *A+  (b- . 5*bA2 ) *A*gam ( 9) + (bA2- . 167*bA3- 
b) *A*gam ( 12 ) ; 

w21=c* ( . 5*bA2-b) *A*gam (7 ) +c*A; 
w22=-l+b*c*A+c* ( .5*bA2-b) *A*gam(8)+c* ( . 167*bA3- 
bA2+b) *A*gam ( 11 ) ; 
w23=0; 

w24=b*c*A+c*A* ( .5*bA2-b) *gam(10) ; 

w25=c*A* ( ( .5*bA2-b) + ( .5*bA2-b) *gam(9) +c* (.167*bA3- 

bA2+b) *A*gam(12) ) ; 

w31= (-bA2+2*b) *A*gam(7) -2*b*A*gam (3) ; 
w32= (-bA2+2*b) *A*gam(8) -2*b*A*gam (5) ; 
w33=-2-2*A-2*b*A*gam(4) ; 

w34=-2*b*A-2*b*A*gam ( 6) +A* (-bA2+2*b) *gam(10) ; 
w35= (-bA2+2*b) *A*gam(9) ; 

w41=b*c*A*gam (3) +d*u ( j ) + (-b*c-b+bA2*c/2 ) *A*gam(7) ; 
w42= (- . 5*bA2+b) *A*gam(ll) + (-b*c- 
b+bA2*c/2) *A*gam (8) +b*c*A*gam (5) ; 
w43=c*A+b*c*A*gam(4) ; 

w44=b*c*A-A-2+b*c*A*gam (6) + (-b*c-b+bA2*c/2 ) *A*gam(10) ; 
w45=-b*A+ (-. 5*bA2+b) *A*gam(12) + (-b*c-b+bA2*c/2) *A*gam(9) ; 
w51=2*c*b*A*gam (7) ; 

w52= (c*bA2-2*c*b) *A*gam(ll) +2*b*c*A*gam (8) +2*d*u ( j ) ; 
w53=0; 

w54=2*b*c*A*gam (10) +2*c*A; 

w55=-2+2*b*c*A+  (bA2-2*b*c) *A*gam(12) +2*b*c*A*gam ( 9) ; 
F=[wll  wl2  wl3  wl4  wl5 
w21  w22  w23  w24  w25 
w31  w32  w33  w34  w35 
w41  w42  w43  w44  w45 
w51  w52  w53  w54  w55] ; 

[E, D] =eig (F) ; 

%eig (F) 

all=exp(D(l,  1) *T)  ; 
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a22=exp(D(2,2) *T) ; 
a33=exp (D (3, 3) *T) ; 
a44=exp (D(4, 4) *T) ; 
a55=exp (D(5, 5) *T) ; 

Pi=[all  a22  a33  a44  a55] ; 

Pl=diag (Pi, 0) ; 

Einv=inv (E) ; 

Ri= [ (all-1) /D (1, 1)  (a22-l)/D(2,2)  (a33-l ) /D (3, 3)  (a44- 

1 ) /D  ( 4 , 4 ) (a55-l)/D(5,5)  ]; 

Rl=diag (Ri, 0) ; 

cl=E*Pl*Einv; 

dl=E*Rl*Einv; 

xl (i+l)=cl (1, l)*xl (i)+cl (1,2) *x2 (i)+cl (1,3) *x3 (i)+cl (1,4) *x 
4 (i) +cl (1, 5) *x5 (i) +d*u (j ) *dl (1,2) ; 

x2 (i+1 ) =cl (2, 1) *xl (i) +cl (2,2) *x2 (i) +cl (2,3) *x3 (i) +cl (2, 4 ) *x 
4 (i) +cl (2, 5) *x5 (i) +d*u ( j ) *dl (2,2) ; 

x3 (i+l)=cl (3, l)*xl (i)+cl (3,2) *x2 (i)+cl (3,3)*x3 (i)+cl  (3, 4) *x 
4(i)+cl(3,5)*x5(i)  +d*u  ( j ) *dl  (3,  2)  ; 

x4 (i+l)=cl (4, 1) *xl (i)+cl (4,2) *x2 (i)+cl (4,3) *x3 (i)+cl (4,4) *x 
4 (i)+cl (4, 5) *x5 (i) +d*u ( j ) *dl (4,2); 

x5 (i+1 ) =cl (5, l)*xl (i)+cl (5,2) *x2 (i)+cl (5,3) *x3 (i)+cl (5, 4) *x 

4 ( i ) +cl (5, 5) *x5 (i) +d*u ( j ) *dl (5,2) ; 

end 

w(j)  = ( (xl  (i+1)  A 3-gam (1 ) *xl  (i+1 ) - 

gam (2) *xl (i+1) A2) /xl (i+1) A3) A2+ ( (xl (i+1) A2*x2 (i+1) - 

gam (3) *xl (i+1) -gam (4) *x2 (i+1) -gam (5) *xl(i+l)A2- 

gam (6) *xl (i+1) *x2 (i+1) ) / (xl (i+1) A2*x2 (i+1) ) ) A2+ ( (xl (i+1) *x2 

(i+1) A2-gam(7) *xl (i+1) -gam (8) *x2 (i+1) - 

gam (9) *xl(i+l)*x2(i+l)- 

gam(10) *x2 (i+1) A2) / (xl (i+l)*x2 (i+1) A2) ) A2+ ( (x2 ( i+1 ) A3- 

gam(ll) *x2 (i+1) -gam (12) *x2 (i+1) A2) /x2 (i+1) A3) A2; 

end 

J=0; 

for  i=l:m 
J=J+w  (i)  ; 
end 
%gama 
J 


function  [t,xl,x2]=discrel6 (T) 
cain=10; 

cas=4 .28 637324003831 6; 
tin=3e02; 

ts=4. 0008 64 742 9034 92e+002; 
v=le-02 ; 
f= . 33e-04 ; 
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z=3e05; 
e=6e04 ; 

R=8 . 314 ; 
h=le04 ; 
rho=le03; 

C=4  ; 

xli= (cain-cas) /cas; 
x2i= (tin-ts) /ts; 
a=v*z/f ; 
b=e/ (R*ts) ; 
c=h*cas/ (rho*C*ts) ; 

qs=-rho*C* (h*v*z*exp ( — b) *cas/ (rho*C) +f* (tin-ts) ) ; 
d=qs/ (f*rho*C*ts) ; 

A=a*exp (-b) ; 

gama(l)=  -0.36378522596406; 

gama (2 ) = 0.09701283041448; 

xl (1) =0; 

x2 ( 1 ) =0; 

x3 (1) =0; 

x4 (1) =0; 

x5  (1 ) =0; 

t (1) =0; 

for  i=l : 10 

j=i-l; 

u= . 9*  (-1) Aj; 

seel=gama (1) *u; 

see2=gama (2) *u; 

gaml (l)=seelA2; 

garni (2) =0; 

gaml (11) =see2A2; 

gaml (12) =0; 

gaml (3) =seel*see2/2; 

gaml (4) =seelA2/2; 

gaml (5) =0; 

gaml ( 6) =0; 

gaml (7) =see2A2/2; 

gaml (8) =seel*see2/2; 

gaml ( 9) =0; 

gaml (10) =0; 

wll=-l-A+ (b- . 5*bA2) *A*gaml (7) ; ; 

wl2=-b*A+ (b- . 5*bA2 ) *A*gaml (8) + (bA2- . 1 67*bA3-b) *A*gaml (11) ; 
wl3=0 ; 

wl4=A* (-b+ (b- . 5*bA2) *gaml (10) ) ; 

wl5=- ( .5*bA2-b) *A+ (b- . 5*bA2 ) *A*gaml (9) + (bA2- . 167*bA3- 
b)  * A* gaml ( 12 ) ; 

w21=c* ( . 5*bA2-b) *A*gaml (7) +c*A; 
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w22=-l+b*c*A+c* ( .5*bA2-b) *A*gaml ( 8 ) +c* ( . 167*bA3- 

bA2+b) *A*gaml (11)  ; 

w23=0; 

w24=b*c*A+c*A* ( . 5*bA2-b) *gaml (10) ; 

w25=c*A*  ( ( .5*bA2-b)  + ( .5*bA2-b) *gaml (9)+c* ( .167*bA3- 
bA2+b) *A*gaml (12) ) ; 

w31= (-bA2+2*b) *A*gaml (7 ) -2*b*A*gaml (3)  ; 
w32= (-bA2+2*b) *A*gaml ( 8 ) -2*b*A*gaml (5) ; 
w33=-2-2*A-2*b*A*gaml (4) ; 

w34=-2*b*A-2*b*A*gaml (6) +A* (-bA2+2*b) *gaml (10) ; 
w35= (-bA2+2*b) *A*gaml (9) ; 

w41=b*c*A*gaml (3) +d*u+ (-b*c-b+bA2*c/2) *A*gaml (7) ; 
w42= (- . 5*bA2+b) *A*gaml (11)+ (-b*c- 
b+bA2*c/2) *A*gaml (8) +b*c*A*gaml (5) ; 
w43=c*A+b*c*A*gaml (4) ; 

w44=b*c*A-A-2+b*c*A*gaml (6) + (-b*c-b+bA2*c/2 ) *A*gaml (10) ; 
w45=-b*A+ (- . 5*bA2+b) *A*gaml (12) + (-b*c-b+bA2*c/2 ) *A*gaml (9) ; 
w51=2*c*b*A*gaml (7) ; 

w52= (c*bA2-2*c*b) *A*gaml (11) +2*b*c*A*gaml (8) +2*d*u; 
w53=0; 

w54=2*b*c*A*gaml (10) +2*c*A; 

w55=-2+2*b*c*A+ (bA2-2*b*c) *A+gaml ( 12 ) +2*b*c*A*gaml (9) ; 
F=[wll  wl2  wl3  wl4  wl5 

w21  w22  w23  w24  w25 

w31  w32  w33  w34  w35 

w41  w42  w43  w44  w45 

w51  w52  w53  w54  w55] ; 

[E,  D] =eig ( F) ; 

%eig (F) 

all=exp (D(l, 1) *T) ; 
a22=exp (D (2, 2) *T) ; 
a33=exp (D (3, 3) *T) ; 
a44=exp (D (4, 4) *T) ; 
a55=exp (D (5, 5) *T) ; 

Pi= [all  a22  a33  a44  a55] ; 

Pl=diag  (Pi, 0) ; 

Einv=inv (E) ; 

Ri=[ (all-1) /D(l, 1)  (a22-l)/D(2,2)  (a33-l) /D (3, 3)  (a44- 

1)/D  (4,4)  (a55-l)/D(5,5)  ] ; 

Rl=diag (Ri, 0) ; 
cl=E*Pl*Einv; 
dl=E*Rl*Einv; 
t (i+l)=i; 

xl (i+1 ) =cl  (1, 1) *xl (i)+cl (1,2) *x2 ( i ) +cl (1,3) *x3 (i) +cl (1, 4) *x 
4 (i)+cl (1, 5) *x5 ( i ) +d*u*dl (1,2); 

x2 (i+1) =cl (2, 1) *xl (i) +cl (2,2) *x2 (i)+cl (2, 3) *x3 (i)+cl (2, 4) *x 
4 (i)+cl (2, 5) *x5 (i) +d*u*dl (2,2) ; 
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x3 (i+l)=cl (3, 1) *xl (i)+cl (3,2) *x2 (i)+cl (3, 3) *x3 (i)+cl (3, 4) *x 
4 (i)+cl (3, 5) *x5 (i) +d*u*dl (3,2); 

x4 (i+l)=cl (4, 1) *xl (i)+cl (4,2) *x2 (i)+cl (4, 3) *x3 (i)+cl (4, 4) *x 
4 (i)+cl (4, 5) *x5 (i) +d*u*dl (4,2) ; 

x5 (i+1) =cl (5, 1) *xl (i)+cl (5,2) *x2 (i)+cl (5, 3) *x3 (i)+cl  (5, 4) *x 

4 (i) +cl (5, 5) *x5 (i) +d*u*dl (5,2) ; 

end 


global  x J 
Jl=1000; 

gama=zeros (4,1)  ; 
options ( 1 ) =0 ; 
options (2) =le-10; 
options (3) =le-10; 
options ( 14 ) =20000; 

optgama=fmins ( ' discn6c ' , gama, options ) 

if ( J<J1)  J1=J; 

Min=Jl ; 
end 

Min 


function  J=discn6c (gama) 
global  x J 
T= . 5 ; 
cain=10; 

cas=4 .286373240038316; 
tin=3e02; 

t s=4. 000864742903492 e+ 002; 

v=le-02; 

f= . 33e-04 ; 

z=3e05; 

e=6e04 ; 

R=8 . 314 ; 
h=le04 ; 
rho=le03; 

C=4 ; 

xli= (cain-cas) /cas; 
x2i= (tin-ts) /ts; 
a=v*z/f ; 
b=e/ (R*ts ) ; 
c=h*cas/ (rho*C*ts) ; 
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qs=-rho*C* (h*v*z*exp (-b) *cas/ (rho*C) +f* (tin-ts) ) ; 
d=qs/ (f*rho*C*ts) ; 

A=a*exp (-b) ; 

u= [ 0 . 9 .8  .7  .6  .5  .4  .3  .2  .1  .0001  -.1  -.2  -.3  -.4  -.5 
.6  -.7  -.8  -0.9] ; 

[n,m]=size (u) ; 
for  j=l:m 
xl (1) =0; 
x2 (1) =0; 
x3 (1) =0; 
x4 ( 1 ) =0; 
x5 (1) =0; 
for  i=l:l 

seel=gama (1) *u ( j ) +gama (2) *u (j ) A2; 

see2=gama (3) *u ( j ) +gama (4) *u ( j ) A2; 

gam(l)=seelA2; 

gam (2 ) =0; 

gam(ll)=see2A2; 

gam (12) =0 ; 

gam (3) =seel*see2/2; 

gam (4) =seelA2/2; 

gam (5) =0; 

gam ( 6) =0; 

gam (7) =see2A2/2; 

gam (8) =seel*see2/2 ; 

gam ( 9) =0; 

gam (10) =0; 

wll=-l-A+ (b- . 5*bA2 ) * A* gam (7) ; ; 

wl2=-b*A+ (b-.5*bA2) *A*gam(8)+ (bA2- . 167*bA3-b) *A*gam(ll) ; 
wl3=0; 

wl4=A* (-b+ (b- . 5*bA2) *gam(10) ) ; 

wl5=- ( .5*bA2-b) *A+ (b- . 5*bA2 ) *A*gam ( 9) + (bA2- . 167*bA3- 
b) *A*gam (12) ; 

w21=c* ( . 5*bA2-b) *A*gam(7) +c*A; 

w22=-l+b*c*A+c* ( .5*bA2-b) *A*gam(8) +c* ( . 167*bA3- 

bA2+b) *A*gam(ll) ; 

w23=0; 

w24=b*c*A+c*A* ( . 5*bA2-b) *gam(10)  ; 

w25=c*A* ( (.5*bA2-b)+ (.5*bA2-b) *gam(9)+c* (.167*bA3- 

bA2+b) *A*gam ( 12 ) ) ; 

w31= (-bA2+2*b) *A*gam (7 ) -2*b*A*gam ( 3 ) ; 
w32= (-bA2+2*b) *A*gam (8) -2*b*A*gam ( 5) ; 
w33=-2-2*A-2*b*A*gam(4) ; 

w34=-2*b*A-2*b*A*gam (6) +A* (-bA2+2*b) *gam(10) ; 
w35= (-bA2+2*b) *A*gam ( 9) ; 

w41=b*c*A*gam (3) +d*u ( j ) + (-b*c-b+bA2*c/2) *A*gam(7) ; 
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w42= (- . 5*bA2+b) *A*gam(ll)  + (-b*c- 
b+bA2*c/2 ) *A*gam ( 8 ) +b*c*A*gam (5); 
w43=c*A+b*c*A*gam(4) ; 

w44=b*c*A-A-2+b*c*A*gam ( 6) + (-b*c-b+bA2*c/2 ) *A*gam(10) ; 
w45=-b*A+ (- . 5*bA2+b) *A*gam(12) + (-b*c-b+bA2*c/2 ) *A*gam ( 9) ; 
w51=2*c*b*A*gam (7) ; 

w52= (c*bA2-2*c*b) *A*gam(ll) +2*b*c*A*gam ( 8 ) +2*d*u ( j ) ; 
w53=0; 

w54=2*b*c*A*gam (10) +2*c*A; 

w55=-2+2*b*c*A+ (bA2-2*b*c) *A*gam ( 12 ) +2*b*c*A*gam ( 9) ; 

F= [wl 1 wl2  wl3  wl4  wl5 

w21  w22  w23  w24  w25 

w31  w32  w33  w34  w35 

w41  w42  w43  w44  w45 

w51  w52  w53  w54  w55] ; 

[E, D] =eig ( F) ; 

%eig (F) 

all=exp(D(l, 1) *T) ; 
a22=exp (D (2, 2) *T)  ; 
a33=exp (D (3, 3) *T) ; 
a44=exp (D (4, 4) *T) ; 
a55=exp (D (5, 5) *T) ; 

Pi= [all  a22  a33  a44  a55] ; 

Pl=diag (Pi,  0) ; 

Einv=inv (E) ; 

Ri=[ (all-1) /D(l, 1)  (a22-l) /D (2, 2)  (a33-l) /D (3, 3)  (a44- 

1 ) /D  ( 4 , 4 ) (a55-l)/D(5,5)  ] ; 

Rl=diag (Ri, 0) ; 

cl=E*Pl*Einv; 

dl=E*Rl*Einv; 

xl  (i+1) =cl  (1, 1) *xl (i)+cl (1,2) *x2 (i)+cl (1, 3)  *x3 (i) +cl (1, 4) *x 
4 (i)+cl  (1, 5)  *x5(i)+d*u(j)  *dl  (1,2)  ; 

x2 (i+1 ) =cl (2, 1) *xl  (i) +cl (2,2) *x2 (i) +cl (2,3) *x3 (i) +cl (2, 4) *x 
4 (i)+cl  (2,5)  *x5  (i ) +d*u  ( j ) *dl  (2,2)  ; 

x3 (i+l)=cl (3,1) *xl (i)+cl (3,2) *x2 (i)+cl (3,3) *x3 (i)+cl (3,4) *x 

4 (i) +cl (3, 5) *x5 (i)+d*u (j) *dl (3,2) ; 

x4 (i+l)=cl (4, 1) +xl (i)+cl (4,2) *x2 (i)+cl (4,3) *x3 (i)+cl (4, 4) *x 

4 (i) +cl (4, 5) *x5 (i) +d*u ( j ) *dl (4,2) ; 

x5 (i+1 ) =cl (5, 1) *xl (i)+cl (5,2) *x2 (i)+cl (5,  3) *x3 (i) +cl (5, 4) *x 

4 (i)+cl (5, 5) *x5 (i) +d*u ( j ) *dl (5,2) ; 

end 

w(j  ) = ( (xl  (i+1)  A3-gam(l)  *xl  (i+1)  - 

gam (2) *xl (i+1) A2) /xl (i+1) A3) A2+ ( (xl (i+1) A2*x2 (i+1 ) - 

gam (3) *xl (i+1) -gam (4) *x2 (i+1) -gam (5) *xl (i+1) A2- 

gam ( 6 ) *xl (i+1) *x2 (i+1) )/ (xl (i+1) A2*x2 (i+1) ) ) A2+ ( (xl (i+1) *x2 

(i+1) A2-gam (7) *xl (i+1) -gam (8) *x2 (i+1 ) — 

gam (9) *xl (i+1 ) *x2 (i+1 ) - 
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gam (10) *x2 (i+1) A2) / (xl (i+1) *x2 (i+1) A2) ) A2+ ( (x2 (i+l)A3- 

gam(ll) *x2 (i+1) -gam (12) *x2 (i+1 ) A2 ) /x2 ( i+1 ) A3) A2; 

end 

J=0; 

for  i=l:m 
J=J+w (i) ; 
end 
%gama 
J 


function  [t,xl,x2]=discrel7 (T) 
cain=10; 

cas=4 .286373240038316; 
tin=3e02 ; 

t s=4. 000864742903492 e+ 002; 

v=le-02; 

f= . 33e-04 ; 

z=3e05; 

e=6e04 ; 

R=8 . 314 ; 
h= 1 e 0 4 ; 
rho=le03; 

C=4 ; 

xli= (cain-cas) /cas; 
x2i= (tin-ts) /ts; 
a=v*  z/f ; 
b=e/ (R* ts ) ; 
c=h*cas/ (rho*C*ts) ; 

qs=-rho*C* (h*v*z*exp (-b) *cas/ (rho*C) +f* (tin-ts) ) ; 
d=qs/ (f*rho*C*ts) ; 

A=a*exp (-b) ; 

gama ( 1 ) = -0.41101893894305; 

gama (2 ) = -0.13007272168497; 

gama (3) = 0.09542670105324; 

gama (4 ) = 0.00330693904224; 

xl ( 1 ) =0; 

x2 ( 1 ) =0 ; 

x3 (1) =0; 

x4 ( 1 ) =0 ; 

x5(l)=0; 

t (1) =0; 

for  i=l:10 

j=i-l; 

u= . 9*  (-1) Aj; 

seel=gama (1) *u+gama (2) *uA2; 
see2=gama (3) *u+gama (4) *uA2; 
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garni  (l)=seelA2; 

gaml (2) =0; 

gaml (ll)=see2A2; 

gaml (12) =0; 

gaml (3) =seel*see2/2; 

gaml (4) =seelA2/2; 

gaml (5) =0; 

gaml ( 6) =0; 

gaml (7) =see2A2/2; 

gaml (8) =seel*see2/2; 

gaml (9) =0; 

gaml (10) =0; 

wll=-l-A+ (b-.5*bA2) *A*gaml (7) ; ; 

wl2=-b*A+ (b- . 5*bA2) *A*gaml ( 8 ) + (bA2- . 1 67*bA3-b) *A*gaml (11) ; 
wl3=0; 

wl4=A*  (-b+ (b- . 5*bA2) *gaml (10) ) ; 

wl5=- ( .5*bA2-b) *A+ (b- . 5*bA2 ) *A*gaml (9) + (bA2-. 167*bA3- 
b)  * A*  gaml  (12)  ; 

w21=c* ( .5*bA2-b) * A* gaml (7) +c*A; 

w22=-l+b*c*A+c* ( .5*bA2-b) *A*gaml (8) +c* ( .167*bA3- 

bA2+b) * A* gaml (11) ; 

w23=0; 

w24=b*c*A+c*A* ( . 5*bA2-b) *gaml (10) ; 

w25=c*A* ( ( .5*bA2-b)+ ( .5*bA2-b) *gaml (9)+c* (.167*bA3- 
bA2+b) *A* gaml (12) ) ; 

w31= (-bA2+2*b) *A*gaml (7) -2*b*A*gaml (3) ; 
w32= (-bA2+2*b) *A*gaml (8) -2*b*A*gaml (5) ; 
w33=-2-2*A-2*b*A*gaml (4) ; 

w34=-2*b*A-2*b*A*gaml (6) +A* (-bA2+2*b) *gaml (10) ; 
w35= (-bA2+2*b) *A*gaml (9) ; 

w41=b*c*A*gaml (3) +d*u+ (-b*c-b+bA2*c/2 ) *A*gaml (7) ; 
w42= (- . 5*bA2+b) *A*gaml (11)+ (-b*c- 
b+bA2*c/2) *A*gaml (8) +b*c*A*gaml (5) ; 
w43=c*A+b*c*A*gaml (4) ; 

w44=b*c*A-A-2+b*c*A*gaml (6) + (-b*c-b+bA2*c/2 ) *A*gaml (10) ; 
w45=-b*A+ (- . 5*bA2+b) *A*gaml (12)  + (-b*c-b+bA2*c/2 ) *A*gaml  (9) ; 
w51=2*c*b*A*gaml (7) ; 

w52= (c*bA2-2*c*b) *A*gaml ( 11 ) +2*b*c*A*gaml (8)+2*d*u; 
w53=0; 

w54=2*b*c*A*gaml (10) +2*c*A; 

w55=-2+2*b*c*A+ (bA2-2*b*c) *A*gaml ( 12 ) +2*b*c*A*gaml (9) ; 
F=[wll  wl2  wl3  wl4  wl5 
w21  w22  w23  w24  w25 
w31  w32  w33  w34  w35 
w41  w42  w43  w44  w45 
w51  w52  w53  w54  w55] ; 

[E,D]=eig  (F)  ; 
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%eig (F) 

all=exp(D(l, 1) *T)  ; 
a22=exp(D(2,2) *T)  ; 
a33=exp (D (3, 3) *T) ; 
a44=exp (D (4, 4) *T) ; 
a55=exp (0(5,5)*!) ; 

Pi= [all  a22  a33  a44  a55] ; 

Pl=diag (Pi, 0) ; 

Einv=inv (E) ; 

Ri=[ (all-1) /D(l, 1)  (a22-l ) /D (2, 2)  (a33-l) /D (3, 3)  (a44- 

1 ) /D ( 4 , 4 ) (a55-l)/D(5,5) ] ; 

Rl=diag (Ri, 0) ; 
cl=E*Pl*Einv; 
dl=E*Rl*Einv; 
t (i+1) =i; 

xl (i+l)=cl (1, 1) *xl (i)+cl (1,2) *x2 (i)+cl (1, 3) *x3 (i)+cl  (1,4) *x 
4 (i)+cl (1, 5) *x5 (i) +d*u*dl (1,2) ; 

x2 (i+1) =cl (2, 1) *xl (i)+cl (2,2) *x2 (i)+cl (2, 3) *x3 (i)+cl (2, 4) *x 
4 (i)+cl (2,5) *x5 (i) +d*u*dl (2,2) ; 

x3 (i+1) =cl (3, 1) *xl (i)+cl (3,2) *x2 (i)+cl (3, 3) *x3 (i)+cl (3, 4) *x 
4 (i)+cl (3, 5) *x5 (i) +d*u*dl (3,2); 

x4 (i+l)=cl (4, 1) *xl (i)+cl (4,2) *x2 (i)+cl (4,3) *x3 (i)+cl (4, 4) *x 
4 (i) +cl (4, 5) *x5 (i) +d*u*dl (4,2) ; 

x5 (i+1 ) =cl (5, 1) *xl (i)+cl (5,2) *x2 (i)+cl (5,3) *x3 (i)+cl (5, 4) *x 

4 ( i ) +cl (5, 5) *x5 (i ) +d*u*dl (5,2) ; 

end 


APPENDIX  D 

MATLAB  CODES  FOR  EXAMPLE  ONE  CHAPTER  3 
(SECOND,  FOURTH  AND  SIXTH  ORDER) 


% CARLMa2.M 

% 


Program_name  = 'CARLMa2.M'  ; 
Revision  = 'Rev.  1.0*  ; 


global  JJ  T w 
format  compact 
Jl=le08; 

[ 0.18733400000000 

-1.51261200000000 
2.38863100000000 
-2.200944] ; 


%coef f 0 
% 

% 

% 


coef f 0 


= randn (4,1) ; 


options  (1) 
%options (2) 
%options (3) 
options (14) 


= 0; 

= le-06; 
= le-06; 
= 5000; 


optcoeff  = fmins (' carlopa2 coef fO, options) ; 


if ( JJ<J1  ) J1  = JJ; 
Min=Jl ; 

end 


% Report  Results 

disp(' Results ') 

disp  ( ' T =' ) , disp (T) 
disp  ( ' wl='),  disp(w(l)) 
dispC  wl2= ' ) , disp  (w  ( 1 ) A2 ) 
disp ( ' w2= ' ) , disp (w (2 ) ) 
disp('  Min  ='),  disp (Min) 
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%disp('  see  = '),  disp(see) 

%disp('  see2  =' ) , disp (see2) 
disp ( ' ' ) 

disp('  Optcoeff  Guesses  '), 

disp([  optcoeff  coeffO  ]) 

disp('  ' ) 

%msg  = sprintf ( ' x tol  = %g',  options (2)  ) ; disp(  msg 

) 

%msg  = sprintf ('  F tol  = %g',  options (3)  ) ; disp(  msg 

) 

msg  = sprintf ('  Max  iter  = %g',  options (14)  );  disp(  msg  ) 

disp(['  File  Name  = ' Program_name  ' ' Revision]) 
disp  ( [ ' Date  = ' date  ] ) 

disp(' End  of  Results ') 


function  JJ=carlopa4 (coeff ) 

global  JJ  see  see2  T w 

T = 15; 
a = coeff (1) ; 
b = coeff (2) ; 
c = coeff (3) ; 
d = coeff (4 ) ; 

maxl=  0; 
minl=-l ; 
max2=  0.5; 
min2=-0 . 5; 
aa  = 0.1; 

H = 5; 
xO  =1; 

cl  = max2-min2; 
dl  = min2 ; 

y = c/ (1+exp (-a) ) +d; 
u =1/ (1  + exp (-y) ) ; 

ul  = cl*u* (1-u) ; 

u2  = cl* (2*uA3-3*uA2+u) ; 

All  = exp (-a) ; 
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A22 

— 

(1+A11)  ; 

yi 

= 

c*b*Al 1/A22 A2 ; 

y2 

= 

c*bA2*Al  1* (All-1) /A22A3; 

fO 

= 

aa*xO+cl*u+dl ; 

al 

= 

aa+ul*yl; 

a2 

= 

u2*yl A2+ul*y2 ; 

A 

= 

[al  a2/2 

2* f 0 2*al ] ; 

bO 

= 

[fO  0] 

wO 

= 

[0  0]  ' ; 

M 

= 

expm (A*T) ; 

I 

= 

eye  (2,2); 

w 

= 

M*wO- ( I-M) *inv (A) *b0; 

cO 

= 

0 . 5*H* [2*x0  1]; 

j 

= 

0 . 5*H*xOA2 ; 

Al 

= 

exp (-y) ; 

A2 

= 

(1+A1)  ; 

A3 

= 

2*cl*Al; 

BO 

= 

0.25* (cl/A2+dl) A2; 

bl 

= 

A3* (cl+dl+dl*Al) /A2A3; 

b2 

= 

A3* (dl*AlA2+A3-cl-dl) /A2A 

Bl 

= 

bl*yl; 

B2 

= 

b2*ylA2+bl*y2; 

B 

= 

0.25* [Bl  B2/2 ] ; 

J1 

= 

cO*w; 

J2 

= 

B0*T; 

J3 

= 

B* (M- I ) * inv (A) *w0; 

J4 

= 

-B* (A*T+I-M) *inv (A) *inv (A) 

J5 

= 

j; 

JJ  = J1+J2+J3+J4+J5; 


% Report  results  to  the  screen 


%clc 
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%disp ( ' 

%disp ( ' J1 ' ) , 
%disp ( ' J2 ' ) / 
%disp ( ' J3 ' ) , 
%disp ( ' J4  ' ) , 
%disp ( ' J5 ' ) , 
disp ( ' JJ' ) , 
%disp ( ' see' ) , 
%disp ( ' see2 ' ) , 
%disp ( ' coeff ' ) , 
disp ( ' ' ) 


disp ( J1 ) 
disp ( J2) 
disp ( J3) 
disp ( J4 ) 
disp ( J5) 
disp ( JJ) 
disp (see) 
disp  (see2 ) 
disp (coeff) 


% CARLMa4.M 

% 


Program_name  = 'CARLMa4.M'  ; 
Revision  = 'Rev.  1.0'  ; 


global  JJ  T w 
format  compact 


Jl=le08; 


%coef f 0 = [0000]'; 

coeff 0 = rand (4,1); 


options (1) 
%options (2 ) 
%options (3) 
options ( 14 ) 


= 0; 

= le-06; 
= le-06; 
= 5000; 


optcoeff  = fmins ( 'carlopa4' , coeff 0, options ) ; 

if ( JJ<J1  ) J1  = JJ; 

Min= J1 ; 
end 


% Report  Results 

disp(' Results ') 

disp  ( ' T = ' ) , disp (T) 
disp('  wl='),  disp(w(l)) 
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disp('  wl2='),  disp (w(l) ^2) 
disp ( ' w2=')/  disp(w(2)) 
disp('  Min  =')/  disp (Min) 

%disp('  see  = ' ) / disp (see) 

%disp('  see2  =' ) , disp (see2) 
disp('  ') 

disp('  Optcoeff  Guesses  '), 

disp([  optcoeff  coeffO  ]) 

disp('  ' ) 

%msg  = sprintf ( ' x tol  = %g',  options (2)  ) ; disp(  msg 

) 

%msg  = sprintf ('  F tol  = %g',  options (3)  ) ; disp(  msg 

) 

msg  = sprintf ('  Max  iter  = %g',  options (14)  );  disp(  msg  ) 

disp(['  File  Name  = ' Program_name  ' ' Revision]) 
disp(['  Date  = ' date  ] ) 

dispC End  of  Results ') 


function  JJ=carlopa4 (coef f ) 

global  JJ  see  see2  T w 

T = 15; 
a = coeff (1) ; 
b = coeff (2) ; 
c = coeff (3) ; 
d = coeff (4 ) ; 

maxl=  0; 
minl=-l ; 
max2=  0.5; 
min2=-0 . 5; 
aa  = 0.1; 

H = 5; 
xO  = 1; 

cl  = max2-min2; 
dl  = min2; 

y = c/ (1+exp (-a) ) +d; 
u =1/ (1+exp (-y) ) ; 
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ul  = cl*u* (1-u) ; 

u2  = cl* (2*uA3-3*uA2+u) ; 

u3  = cl* (-6*uA4+12*uA3-7*uA2+u) ; 

u4  = cl* (24*uA5-60*uA4+50*uA3-15*uA2+u) ; 

All  = exp (-a) ; 

A22  = (1+A11); 

yl  = c*b*All/A22A2 ; 

y2  = c*bA2*All* (All-1) /A22A3; 

y3  = c*bA3*All* (All A2-4*A11+1 ) /A22A4; 

y4  = c*bA4*All* (All-1) * (A11A2-10*A11+1 ) /A22A5 

fO  = aa*xO+cl*u+dl ; 
al  = aa+ul*yl; 
a2  = u2*ylA2+ul*y2; 
a3  = u3*ylA3+3*u2*yl*y2+ul*y3; 

a4  = u4*ylA4+6*u3*y2*yl A2+4*u2*yl*y3  . . . 

+3*u2*y2A2+ul*y4; 

A = [al  a2/2  a3/6  a4/24 

2*f 0 2*al  a2  a3/3 
0 3* f 0 3*al  3*a2/2 

0 0 4*f 0 4*al  ] ; 

bO  = [fO  0 0 0] ' ; 
wO  = [0000]'; 

M = expm(A*T); 

I = eye  (4,4); 
w = M*wO- (I-M) *inv(A) *bO; 
cO  = 0 . 5*H* [2*x0  100]; 
j = 0 . 5*H*xOA2 ; 

Al  = exp(-y); 

A2  = (1+A1) ; 

A3  = 2*cl*Al ; 

BO  = 0.25* (cl/A2+dl) A2; 

bl  = A3* (cl+dl+dl*Al) /A2A3; 

b2  = A3* (dl*Al A2+A3-cl-dl ) /A2A4 ; 

b3  = A3* (dl*AlA3+ (4*cl-3*dl) *A1A2+ (-7*cl- 

3*dl ) *Al+cl+dl ) /A2A5; 

b4  = A3* (dl*AlA4+ (8*cl-10*dl) *A1 A3-33*cl*Al A2 

+ (18*cl  + 10*dl) *Al-cl-dl)  /A2A6; 

Bl  = bl*yl; 
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B2  = b2*ylA2+bl*y2; 

B3  = b3*ylA3+3*b2*yl*y2+bl*y3; 

B4  = b4*yl A4+6*b3*y2*ylA2+4*b2*yl*y3  ... 
+3*b2*y2A2+bl*y4 ; 

B = 0.25* [B1  B2/2  B3/6  B4/24]; 

J1  = c0*w; 

J2  = B0*T; 

J3  = B* (M-I ) *inv (A) *w0 ; 

J4  = -B* (A*T+I-M) *inv (A) *inv (A) *b0; 

J5  = j ; 

JJ  = J1+J2+J3+J4+J5; 

see  = abs ( (a2*w (2 ) +a3*w (3) +a4*w (4 ) ... 

-a2*w(l) A2-a3*w(l) A3-a4*w(l) A4) / . . . 

( f 0+al*w ( 1 ) +a2*w ( 1 ) A2+a3*w ( 1 ) A3+a4*w ( 1 ) A4 ) ) ; 
see2=  abs ( (w (2 ) - 

w(l)A2)/w(l)A2) +abs ( (a2*w (2) +a3*w (3) +a4*w (4 ) ... 

-a2*w(l) A2-a3*w(l) A3-a4*w(l) A4) / . . . 

(a2*w ( 1 ) A2+a3*w ( 1 ) A3+a4*w ( 1 ) A4 ) ) ; 


%wl=w(l) 

% Report  results  to  the  screen 


%clc 

%disp ( ' 

%disp ( ' J1 ' ) , 
%disp ( ' J2 ' ) / 
%disp ( ' J3 ' ) , 
%disp ( ' J4 1 ) , 
%disp ( ' J5 ' ) , 
disp ( ' JJ' ) , 
%disp ( ' see ' ) , 
%disp ( ' see2 ' ) , 
%disp ( ' coef f ' ) , 
disp  ( ' ' ) 


disp ( J1 ) 
disp ( J2) 
disp ( J3) 
disp ( J4) 
disp ( J5) 
disp ( JJ) 
disp (see) 
disp (see2) 
disp (coef f ) 


% CARLMa  6 . M 

% 

Program  name  = 


'CARLMa6.M'  ; 
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Revision  = 'Rev.  1.0' 

global  JJ  T w 
format  compact 


Jl=le08 ; 

%coef f 0 = [0000]'; 

coeffO  = randn(4,l); 


options ( 1 ) 
%options (2) 
%options (3) 
options (14) 


= 0; 

= le-06; 
= le-06; 
= 5000; 


optcoeff  = fmins (' carlopa6 ', coeffO, options) ; 


if ( JJ<J1  ) J1  = JJ; 
Min= J1 ; 

end 


% Report  Results 

disp(' Results ') 

disp('  T ='),  disp(T) 
disp ( ' wl='),  disp(w(l)) 
disp('  wl2='),  disp(w(l)A2) 
disp('  w2='),  disp(w(2)) 
disp ( ' Min  = '),  disp (Min) 

%disp('  see  = '),  disp (see) 

%disp('  see2  =' ) , disp (see2) 
disp('  ') 

disp('  Optcoeff  Guesses  ') 

disp([  optcoeff  coeffO  ]) 

disp ( ' ' ) 

%msg  = sprintf ( ' x tol  = %g',  options (2)  ) 
) 

%msg  = sprintf ('  F tol  = %g',  options (3)  ) 

) 

msg  = sprintf ('  Max  iter  = %g',  options (14)  ) 

disp(['  File  Name  = ' Program_name 
disp(['  Date  = ' date  ] ) 


; disp ( msg 
; disp(  msg 
; disp(  msg  ) 


I V 


Revision] ) 
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disp(' End  of  Results ') 


function  JJ=carlopa6 (coef f ) 

global  JJ  T w 

T = 15; 
a = coef f ( 1 ) ; 
b = coef f (2 ) ; 
c = coeff (3) ; 
d = coeff (4) ; 

maxl=  0; 
minl=-l ; 
max2=  0.5; 
min2=-0 . 5; 
aa  =0.1; 

H = 5; 
xO  = 1; 

cl  = max2-min2 ; 
dl  = min2 ; 

y = c/ (1+exp (-a) ) +d; 
u =1/ (1+exp (-y) ) ; 

ul  = cl*u* ( 1-u) ; 
u2  = cl* (2*uA3-3*uA2+u) ; 
u3  = cl* (-6*uA4+12*uA3-7*uA2+u) ; 
u4  = cl* (24*uA5-60*uA4+50*uA3-15*uA2+u) ; 
u5  = cl* (-120*uA6+360*uA5-390*uA4+180*uA3-31*uA2+u) ; 
u6  = cl* (720*uA7-2520*uA6+3360*uA5-2100*uA4+602*uA3- 
63*uA2+u) ; 

All  = exp (-a) ; 

A22  = ( 1+A11 ) ; 

yl  = c*b*All/A22A2; 
y2  = c*bA2*All* (All-1) /A22A3; 
y3  = c*bA3*All* (A11A2-4*A11+1) /A22A4; 
y4  = c*bA4*All* (All-1) * (A11A2-10*A11+1) /A22A5; 
y5  = c*bA5*All* (A11A4-26*A11A3+66*A11A2-26*A11+1) /A22A6; 
y6  = c*bA6*All* (All-1) * (A11A4-56*A11A3+246*A11A2- 
5 6 * A1 1 + 1 ) / A2  2 A 7 ; 


104 


fO  = aa*xO+cl*u+dl ; 
al  = aa+ul*yl; 
a2  = u2*ylA2+ul*y2; 
a3  = u3*ylA3+3*u2*yl*y2+ul*y3; 

a4  = u4*ylA4+6*u3*y2*ylA2+4*u2*yl*y3  . . . 

+3*u2*y2A2+ul*y4 ; 

a5  = u5*ylA5+10*u4*y2*ylA3+10*u3*y3*ylA2  ... 

+15*u3*yl*y2A2+10*u2*y2*y3+5*u2*yl*y4+ul*y5; 
a6  = u6*yl A6+15*u5*y2*yl A4+20*u4*ylA3*y3  ... 
+45*u4*ylA2*y2A2+15*u3*ylA2*y4  . . . 
+60*u3*yl*y2*y3+15*u3*y2A3+10*u2*y3A2  . . . 
+15*u2*y2*y4+6*u2*yl*y5+ul*y6; 


al 

a2/2 

a3/6 

a4/24 

a5/120 

a6/720 

2*  f 0 

2*al 

a2 

a3/3 

a4/12 

a5/60 

0 

3*  f 0 

3*al 

3*a2/2 

a3/2 

a4/8 

0 

0 

4*  f 0 

4*al 

2*a2 

2*a3/3 

0 

0 

0 

5*  f 0 

5*al 

5*a2/2 

0 

0 

0 

0 

6*  f 0 

6*al 

bO  = [fO  0 0 0 0 0] ' ; 
wO  = [000000]'; 

M = expm(A*T); 

I = eye (6, 6) ; 
w = M*wO- (I-M) *inv (A) *bO; 
cO  = 0 . 5*H* [2*x0  10000]; 
j = 0 . 5*H*xOA2 ; 

Al  = exp ( -y ) ; 

A2  = (1+A1) ; 

A3  = 2*cl*Al; 

BO  = 0.25* (cl/A2+dl) A2; 

bl  = A3* (cl+dl+dl*Al) /A2A3; 

b2  = A3* (dl*AlA2+A3-cl-dl) /A2A4; 

b3  = A3* (dl*AlA3+ (4*cl-3*dl) *A1A2+ (-7*cl- 

3*dl) *Al+cl+dl) /A2A5; 

b4  = A3* (dl*AlA4+ (8*cl-10*dl) *AlA3-33*cl*Al A2  ... 

+ (18*cl  + 10*dl) *Al-cl-dl) /A2A6; 
b5  = A3* (dl*AlA5+ ( 16*cl-25*dl ) *A1A4+ (- 

131*cl+40*dl ) *A1A3  ... 

+ (171*cl+40*dl) *A1A2+ (-41*cl-25*dl) *Al+cl+dl) /A2A7; 
b6  = A3* (dl*AlA6+ (32*cl-56*dl) *A1A5+ (- 

473*cl+245*dl) *A1A4  ... 

+1208*cl*AlA3+ (-718*cl- 
245*dl) *A1A2+ (88*cl+56*dl) *Al-cl-dl)  /A2A8; 
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B1  = bl*yl; 

B2  = b2*ylA2+bl*y2 ; 

B3  = b3*ylA3+3*b2*yl*y2+bl*y3; 

B4  = b4*yl A4+6*b3*y2*ylA2+4*b2*yl*y3  ... 
+3*b2*y2A2+bl*y4 ; 

B5  = b5*yl A5+10*b4*y2*ylA3+10*b3*y3*yl A2  ... 

+15*b3*yl*y2A2+10*b2*y2*y3+5*b2*yl*y4+bl*y5; 
B6  = b6*ylA6+15*b5*y2*ylA4+20*b4*yl A3*y3  ... 

+45*b4*yl A2*y2A2+15*b3*ylA2*y4  . . . 
+60*b3*yl*y2*y3+15*b3*y2A3+10*b2*y3A2  . . . 
+15*b2*y2*y4+6*b2*yl*y5+bl*y6; 

B = 0.25* [B1  B2/2  B3/6  B4/24  B5/120  B6/720] ; 

J1  = c0*w; 

J2  = B0*T; 

J3  = B* (M-I ) *inv (A) *w0; 

J4  = -B* (A*T+I-M) *inv (A) *inv (A) *b0; 

J5  = j ; 

JJ  = J1+J2+J3+J4+J5; 


%wl=w(l) 


% Report  results  to  the  screen 


%clc 

%disp ( ' 

%disp ( ' J1 ' ) / 
%disp ( ' J2 ' ) , 
%disp ( ' J3 ’ ) / 
%disp ( ' J4 ' ) , 
%disp ( ' J5 ' ) , 
disp ( ' JJ' ) , 
%disp ( ' see ' ) , 
%disp ( ' see2 ' ) , 
%disp ( ' coef f ' ) , 
disp ( ' ' ) 


disp ( J1 ) 
disp ( J2) 
disp ( J3) 
disp ( J4 ) 
disp ( J5) 
disp ( JJ) 
disp (see) 
disp (see2) 
disp (coef f ) 


ul=[]  ; 
tl=  [ ] ; 
U = [ ] ; 
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J1=0; 
xx=l  ; 


tO  = 0; 
tf  = 15; 
xO  = 0; 
tol=le-15; 

a = [ 0.41167822218513 

0.02250442166144 
-0.62973951397489 
-0.14818303119785] ; 

aa  = 0.1; 
mx2=  0.5; 
mn2=-0 . 5; 


[t,y]  = ode45 ( ' carlint2 ' , tO,  tf , xO,  tol) ; 

[m,n]  = size (y) ; 

for  i=l:m 

if  i==l 

dt=t  (i)  ; 
else 

dt=t (i)-t(i-l) ; 

end 

zl  = a (3) / (1+exp (- (a (1) +a (2) *y (i) ) ) ) +a  (4) 

%zl  = a (3) / (1+exp (- (a (1) +a (2) *y (i) ) ) ) +a  (4) 
% +a (7) / (1+exp (- (a (5) +a (6) *y  (i) ) ) ) ; 

%zl  = a (3) / (1+exp (- (a (1) +a (2) *y (i) ) ) ) +a (4) 
% +a (7) / (1+exp (- (a (5) +a (6) *y (i) ) ) ) ... 

% +a (10) / (1+exp (— (a (8)+a (9) *y(i)  ) ) ) ; 

%zl  = a (3) / (1+exp (- (a (1) +a (2) *y(i) ) ) )+a(4) 
% +a (7) / (1+exp (- (a (5) +a (6) *y (i) ) ) ) ... 

% +a (10) / (1+exp (- (a (8) +a (9) *y (i) ) ) ) .. 

% +a (11) / (1+exp (— (a (12) +a (13) *y(i) ))); 

ul (i)  = (mx2-mn2) / (1+exp (-zl) ) +mn2; 

J1  = J1  + ul (i) A2*dt; 

U = [ U ; ul  (i)  ] ; 


end 
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J2  = 5* (y (m)+xx) A2/2; 

J = J2  + Jl/4; 

subplot (211 ) , plot(t,ul) 
axis ([0  15  -0.25  0]) 
xlabel ( ' t ' ) 
ylabel ('u' ) 

title (' 6th  order  No.  Nodes=  1') 

subplot (212) , plot(t,y) 
xlabel ( ' t ' ) 
ylabel ( 'x' ) 


function  yprime=carlint2 ( t, x) 


a = [ 0.41167822218513 

0.02250442166144 
-0.62973951397489 
-0.14818303119785] ; 


aa  = 0.1; 
mx2=  0.5; 
mn2=-0 . 5; 


z=a (3) / (1+exp (- (a (1) +a (2) *x) ) ) +a (4) ; 

%z=a(3)/ ( 1+exp  (- (a (1 ) +a (2 ) *x) ) )+a (4)+a (7) / (1+exp ( 
(a(5)+a(6)*x)  ) ) ; 

%z=a (3) / (1+exp (- (a (1) +a (2) *x) ) ) +a (4) +a (7) / (1+exp ( 
(a (5) +a (6) *x) ) ) ... 

% +a (10) / (1+exp (- (a (8) +a (9) *x) ) ) ; 

%z=a (3) / (1+exp (— (a (l)+a (2) *x) ) )+a(4)+a(7)/ (1+exp ( 
(a (5) +a (6) *x) ) ) ... 

% +a (10) / (1+exp (-(a(8)+a(9)*x) ) ) +a ( 11 )/( 1+exp (- 
(a (12) +a (13) *x) ) ) ; 

yprime (1) =aa* (x+1) + (mx2-mn2) / ( 1+exp (-z) ) +mn2; 


dt  = 0.05; 
T = 15; 
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N = T/dt; 

Y = [ ] ; 

Th  = [ ] ; 

U = [ ] ; 

X = [ ] ; 

J1  = 0; 

a =0.1; 

H = 5; 

for  i=l:N 

if  i==l 

tO  = 0; 
tf  = dt; 
xO  = 1; 

else 

tO  = ( i— 1 ) *dt; 
tf  = i*dt; 
xO  = y (m) ; 

end 


tol=  le-10; 

[ t, y] =ode45 ( ' exaoptl ' , tO, tf , xO, tol) ; 

[m, n] =size (y) ; 

Y = [ Y ; y (m)  ] ; 

Th=  [ Th;  t (m) ] ; 


yl  = a*  (T-t  (m)  ) ; 

K (m)  = inv (exp (-yl ) -H* (exp (-yl) -exp (yl ) ) /a) *H*exp (yl ) ; 

u (m)  = -2*K (m) *y (m) ; 

Yl  (m)  = y (m)  -1; 

J1  = J1  + u(m)A2*dt; 
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U = [ U ; u (m)  ] ; 

X = [ X ; Y1  (m)  ] ; 

end 

J = 5*y (m) A2/2  + Jl/4; 


subplot  (211 ) , plot(Th,U) 
title ('  Exact  Solution  ') 
%axis ( [0  10  -0.4  -0.1] ) 
xlabel ( ' t ' ) 
ylabel ( 'u' ) 

subplot  (212) , plot(Th,X) 
%axis ( [0  5 -0.08  0] ) 
xlabel  ( ' t ' ) 
ylabel ( ' x ' ) 

%subplot  (212) , plot ( t , K) 
%xlabel  ( ' t ' ) 

%ylabel ( 'K' ) 


function  xprime=exaoptl ( t , x) 

T = 15; 
a = 0.1; 

H = 5; 

yl=  a* (T-t) ; 

K = inv (exp ( — y 1 ) -H* (exp ( — y 1 ) -exp (yl ) ) /a) *H*exp (yl ) ; 
u = -2*K*x; 


xprime (1) =a*x+u; 


APPENDIX  E 

CALCULATIONS  RELATED  TO  LEMMA  1 
(CHAPTER  2) 


Case  E-l : n = 2 and  k = 2 


B = 


b„  bI2 

L^21  ^22 


B®  B = 


2b  n 
b2i 

b21 
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u,2 

b„  +b22 


u12 
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b„  +b22 
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'12 


'12 


2b 


22  . 
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1 0 0 
0 1 0 
0 1 0 
0 0 1 


[b®b]r  = 


bn 

2b12 

0 

21 

bll  +b22 

b12 

21 
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b12 

0 

2b2i 
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Case  E-2: 

b,, 


B = 


'21 
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Case  E-3 : n = 3 and  k = 2 
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